Methods and systems for brain function analysis

ABSTRACT

Various methods and systems are provided for cerebral diagnosis. In one example, among others, a method includes obtaining EEG signals from sensors positioned on a subject; conditioning data from the EEG signals to remove artifacts; generating a cerebral network model based at least in part upon the conditioned data; determining network features based upon the cerebral network model; and determining a cerebral condition of the subject based at least in part upon the network features. In another example, a method includes determining a recording condition of a positioned EEG sensor and providing an indication of an unacceptable recording condition of the EEG sensor. In another example, a system includes an EEG recording module to acquire signals; a signal conditioning module to condition signal data; a signal analysis module to determine signal features and cerebral network features; and a condition classification module to determine a cerebral condition of the subject.

CROSS REFERENCE TO RELATED APPLICATIONS

This application claims priority to co-pending U.S. provisionalapplication entitled “Methods and Systems for Cerebral Diagnosis” havingSer. No. 61/612,647, filed Mar. 19, 2012, the entirety of which ishereby incorporated by reference.

BACKGROUND

The human brain is a very delicate organ that makes possible complexbehavioral decision making. Information processing within the humanbrain is so sophisticated and complex that it cannot be accessedentirely even with the advanced technology available today. Once a brainis damaged, it is often very hard to achieve full recovery. Theimportance of accurate, timely diagnoses of brain abnormality is crucialin many clinical settings including the emergency room (ER) or intensivecare unit (ICU). However, most mental and neurological states areevaluated mainly through interviews and subjective exams based on thesubjects' temporary performance at that time. There is no objectivequantitative test for evaluating baseline brain function. Imagingtechnologies such as standard magnetic resonance imaging (MRI) show onlystructure within the brain without providing an indication of dynamicbrain function. EEG is the most effective method for evaluating brainfunction, but interpretation requires interpretation of multichannelgraphs based on visual analysis by highly trained experts.

BRIEF DESCRIPTION OF THE DRAWINGS

Many aspects of the present disclosure can be better understood withreference to the following drawings. The components in the drawings arenot necessarily to scale, emphasis instead being placed upon clearlyillustrating the principles of the present disclosure. Moreover, in thedrawings, like reference numerals designate corresponding partsthroughout the several views.

FIG. 1 is a block diagram illustrating an example of a system forevaluating a condition of a brain in accordance with h variousembodiments of the present disclosure.

FIGS. 2A and 2B are a flowchart illustrating an example of functionalityof the system of FIG. 1 in accordance with various embodiments of thepresent disclosure.

FIGS. 3A and 3B are a flowchart illustrating examples of featureextraction, network modeling, and classification of the flowchart ofFIG. 2B in accordance with various embodiments of the presentdisclosure.

FIG. 4 is a graphical representation of an example of a processor systemsuitable for implementing the system of FIG. 1 in accordance withvarious embodiments of the present disclosure.

DETAILED DESCRIPTION

Disclosed herein are various embodiments of methods and systems relatedto diagnosis of cerebral conditions which cause disturbances in brainfunction. Reference will now be made in detail to the description of theembodiments as illustrated in the drawings, wherein like referencenumbers indicate like parts throughout the several views.

Electroencephalography is a technology for measuring the voltage andfrequency of electrical activity from neurons in the cerebral cortex.Electroencephalogram (EEG) electrodes can record brainwaves usingelectrodes attached to the scalp or, through electrodes placed on thesurface of the brain (subdural electrodes) or within brain tissue (depthelectrodes) using surgical procedures. A scalp EEG is a non-invasiveprocedure which provides useful information about brain state andfunction. This methodology is used in many fields of neuroscience (e.g.,psychology, epilepsy, brain machine interface, and sleep research) forrecording and analyzing brain state and function. It is used widely as adiagnostic tool in clinical neurology to evaluate and monitor brainfunction and to identify disturbances in the function of the braincaused by a variety of insults to the brain, such as concussion,traumatic injury, stroke, tumor, encephalopathies due to toxins ormetabolic disturbances and seizures. Many disturbances of brain functioncan be identified by analysis of brief multichannel EEG recordings usingelectrodes placed in specific locations on the scalp, based of referenceanatomical landmarks. The most widely accepted system of electrodeplacement is the International 10/20 System of electrode placement. Byanalyzing the multi-channel EEG data of a short period of time, ongoingbackground activity can be assessed. Standard routine EEG recordings areapproximately 20 minutes in duration. However, background EEG can beassessed with only a few seconds of EEG recording, if the state of thesubject is known (e.g., alert, drowsy or light sleep states). A normalbrain generates signals with characteristic frequencies, waveforms andspatial organization. Normal brain electrical activity is remarkablysymmetrical over the two cerebral hemispheres. In clinical practice, theEEG is analyzed visually to detect diffuse bilaterally disturbances inongoing background activity as well as focal or lateralizeddisturbances. This information provides useful diagnostic insight. TheEEG is useful for evaluation of chronic conditions, such as dementia.However, its most important use is in the evaluation of acute orsubacute conditions presenting as altered mental status or impairedsensorium, such as stupor or coma.

Analysis and interpretation of EEG recordings is performed by experts(usually neurologists with training in the interpretation of EEGrecordings) based on visual inspection of multichannel recordingsdisplayed as multichannel graphs of signal voltage over time. Thespatial temporal patterns of brain electrical activity can be analyzedthough quantitative analysis of the spatial and temporal properties ofthis activity, generating a mathematical model of the activity,analyzing the properties of the model, and comparing those properties toa standard of norms based on the properties of normal individuals ofsimilar age. A useful mathematical model for analyzing brain state andfunction is a network model, based on graph theory. Mathematicalanalysis of local and global properties allow identification of normalphysiological states and allow identification of focal, lateralized ordiffuse disturbances in physiological function and states. While briefEEG recordings provide useful diagnostic information regarding functionand state at the time of recording, analysis of long-term recording isuseful for monitoring brain conditions such as level of consciousness,identifying abnormal transients in the signal, such as interictalepileptiform discharges, seizures, acute cerebral ischemia, hypoglycemiaor anoxia. Analyzing the spatial-temporal dynamics of long-term EEGrecording data can be achieved through use of quantitative dynamicalnetwork models. This approach can provide diagnostic informationpertaining to the physiological states of the brain, and can be used toidentify transient pathological conditions. The above described approachto network analysis of brain electrical activity can be achieved throughcomputer-based algorithms designed to identify artifact, condition thesignal, generate quantitative measures of signal properties from each ofmultiple EEG channels derived from multiple electrodes placed instandard locations on the scalp, generating a network model, calculatinglocal and global properties of the network and comparing to standardnetwork norms derived from normal subjects. These network properties canbe monitored over time and compared with the patient or subject'sbaseline to detect significant changes in state or development oftransient pathological conditions. The algorithm can provide detailedquantitative output or can summarize the results and categorize them asnormal or abnormal. Abnormalities can be categorized as to focal,lateralized or diffuse, and the anatomical location focal andlateralized abnormalities can be reported. This output can be written asa report or depicted as a graph. The algorithm can be trained andcompared to the interpretation of expert electroencephalographers tooptimize the sensitivity and specificity of the algorithm. An appendix,which is hereby incorporated by reference in its entirety, providesadditional information about EEG analysis of brain dynamical behavior.

This disclosure presents systems and methods for reliably evaluating arecorded EEG in real time by algorithms and providing an immediateindication of the cerebral condition. The acquired EEG data may beevaluated in real time and/or may be stored for subsequent evaluation.The system allows EEG data to be recorded and evaluated when no EEGtechnical personnel or neurological interpreters are readily available,allowing its use as a screening tool to assist physicians when suchpersonnel are not available. The original EEG data may also be storedfor subsequent visual analysis and/or may be transmitted to remote sitesfor review and interpretation by experts.

EEG data can be recorded with small, portable, inexpensive instrumentsthat do not require special shielded facilities or the subject remainingmotionless for long periods of time. Thus, an EEG can be utilized innoisy point-of-care environments, such as Emergency rooms and IntensiveCare Units and with subjects who may be uncooperative. Portable EEGunits can also be utilized in emergency vehicles and in the field.Referring to FIG. 1, shown is a block diagram illustrating an example ofa system 100 for evaluating a condition of a brain. In the embodiment ofFIG. 1, the system 100 includes an electrode application module 103, anEEG recording module 106, a signal conditioning module 109, a signalanalysis module 112, and a condition classification module 115. Thesystem 100 may contain an interactive display which may provide, e.g.,step-by-step instructions on electrode placement, establishingconnections, preparing the subject and initiating the recording, etc.The system 100 allows for rapid acquisition and analysis of EEG signals,measurement of the spatial-temporal characteristics of the signal,analysis of local, regional, and diffuse signal characteristics,characterization of network features of EEG signals, determination ofwhether or not these EEG characteristics are normal or abnormal, andclassification of abnormal EEG recordings to determine whether theabnormalities are local, lateralized, multifocal or diffuse. The resultsmay be displayed as in a text and/or graphical format that may beavailable for immediate use by emergency room personnel, intensive carepersonnel, and emergency medical technicians.

Referring to FIGS. 2A-2B, shown is a flowchart illustrating variousfunctions that may be implemented by modules of the system 100. As shownin FIG. 2A, after operation of the system has been initiated (e.g., byturning on the device), an electrode application procedure may beprovided in box 203 by the electrode application module 103 forrendering on the system display. The EEG provides direct informationabout brain functions through analysis of brain electrical activity.Depending on the location and the type of the electrodes, EEG signalscan reveal different levels of neuronal activity. Numerous scalp EEGelectrodes may be applied to a subject to obtain EEG data containinginformation from brain activity in both temporal and spatial domains.The EEG electrodes can include individual electrodes and/or an array ofelectrodes that are positioned on the scalp. The location of electrodeplacement on the scalp can follow, e.g., a 10-20 system or otherappropriate system as can be appreciated. In the 10-20 system, distinctlandmarks of the subject's head are first identified and then electrodesare placed at 10% or 20% distance intervals along the landmarks. In someimplementations, intracranial electrodes may be utilized. The recordingintegrity and impedance of the electrodes are checked in box 206 todetermine if there is a problem with the placement and/or operation ofthe EEG electrodes. If there is a problem, the problematic electrodesmay be indicated in box 209 on the system display and the appropriateapplication or correction procedure(s) may be provided in box 203.Guiding the operator through simple set-up and operating procedures toobtain a technically adequate EEG recording reduces evaluation errors.

If the electrode integrity and impedance is acceptable in box 206, thenan EEG is obtained under specified conditions in box 212 by the EEGrecording module 106. For example, EEGs may be obtained with thesubject's eyes open and closed or other specified conditions of thesubject. Instructions may be rendered for display on the system displayto prompt the specified condition (or action) in the monitored subject.For example, during a routine EEG examination, a subject is usuallyasked to relax and open eyes for a period of time and then close eyesfor another period of time. The acquired EEG signals are then processedin box 215. For example, amplification and filtering may be applied toenhance the signal-to-noise ratio (SNR) of the EEG signals. Analog EEGsignals from the electrodes may also be digitized for communication andstorage of the information. In box 218, acceptability of the recordingquality of the EEG data is confirmed. For example, all channels of theprocessed EEG signal may be analyzed for the presence of excessiveartifacts that may contaminate the EEG data. Criteria for acceptablesignal quality may be predefined to ensure acceptable electrode contact,electrode impedance, and minimal contamination by common artifacts.

If the EEG data is not acceptable, then the system can return to box 206to recheck electrode integrity and impedance. Common technical problemsthat degrade the recording (e.g., excess muscle or movement artifacts)may also be determined in box 218. Instructions may be provided throughthe system display to guide the operator in methods to eliminate orattenuate those artifacts before repeating the acquisition of EEGsignals in box 212. Subsequent recorded EEG data may be re-evaluated andthe operator notified of persistent problems, at which time the operatormay attempt to obtain further EEG signals or may abort the procedure. Ifthe EEG data is acceptable, the digitized data can be stored in a datastore or other memory in box 221. The stored EEG data may be transmittedthrough a wireless or wired network connection (e.g., cellular,Bluetooth, Ethernet, etc.) for remote evaluation, analysis, and/orconfirmation.

The acceptable EEG data is further processed and/or filtered by thesignal conditioning module 109 to remove common recording artifacts asillustrated in FIG. 2B. For example, eye movement artifacts may bedetected and removed in box 224, electromyogram (muscle movement)artifacts may be detected and removed in box 227, and electrode relatedartifacts such as, e.g., electromagnetic interference from nearbyinstruments may be detected and removed in box 230. Other artifacts suchas, e.g., 60 Hz line signals and signals produced by mechanicalventilators and other instruments may also be detected and removed fromthe EEG data by the signal conditioning module 109. For example, epochsof data may be examined sequentially for the presence of an artifact. Ifa segment contains an excessive artifact, it may then be excluded fromsubsequent analysis. After the signal conditioning module 109 identifiesand discards the artifact-contaminated segments, the remaining EEG datais evaluated in box 233 to ensure that a useable signal of sufficientduration has been acquired and is available for analysis An initialinterpretation may be generated based upon a brief EEG recording (e.g.,2 to 5 minutes). If the available data is not long enough for reliableanalysis, the system 100 will inform the operator and return to box 212to obtain additional EEG data. If the available EEG data is long enoughfor analysis or evaluation, the EEG data is processed by the signalanalysis module 112 to extract features in box 236.

The system 100 may also provide an option to continue recording toobtain a complete routine EEG (typically around 20-30 minutes ofrecording) or for continued monitoring the EEG for changes in brainfunction, such as intermittent seizures, diffuse or focal ischemia, andchanges in alertness or level of consciousness. In otherimplementations, the system 100 may be placed in a monitoring mode inwhich epochs of the EEG are analyzed as they are acquired to detecttransient abnormalities or state changes in the subject. In themonitoring mode, a continuous or intermittent analysis may be providedgraphically and/or a summary report may be provided intermittently atspecified intervals. For example, the interval between reports can be adefault interval (e.g., every 10 minutes) or can be an interval that isselected by the operator.

The extracted features of the multichannel EEG data from box 236 may beused to provide a quantitative description of the spatiotemporalcharacteristics of the signals, including local and regionalcharacteristics, inter-hemispheric asymmetries, and local and globalnetwork connectivity characteristics. The extracted features may be,e.g., linear, non-linear, univariate, or bivariate statistics. Theextracted features from box 236 may be provided as inputs for networkmodeling in box 239 and for classification of the cerebral condition inbox 245 of the condition classification module 115. For example, if thefeature is univariate, such as entropy, each EEG channel will have afeature time series. The network modeling in box 239 is implementedbased upon the degree of association between these univariate featuresbetween channels. If the feature is bivariate (or a relationship betweentwo channels), the network modeling in box 239 may be directlyimplemented through the values of bivariate features. After the networkmodel has been constructed in box 239, network features may be extractedin box 242 from the network model and provided to the conditionclassification module 115 for classification in box 245. Classificationof the cerebral condition can be more precise by including the extractednetwork features in the evaluation.

Classification of the cerebral condition in box 245 may be based, atleast in part, upon comparison of extracted features from boxes 236 and242 by comparison with established norms to determine if they indicate anormal condition within normal limits or an abnormal condition. Inaddition being able to utilize correlations between specific EEGfindings and pathologies, the condition classification module 115analyzes EEG signals through a network perspective. The functionalnetwork reflects the connectedness among brain regions in terms ofneuron activity. The brain functional network may be represented as agraph by defining vertices and edges in the EEG data. If the EEGchannels are designated as the vertices of a graph, an edge between twovertices signifies a functional connection between two EEG channels. Alarger correlation between two EEG channels indicates the presence of anedge between the channels. Edges may also be values quantifying how wellthe two vertices correlate in weighted graphs. Applying graphtheoretical analysis to EEG data can reveal topological characteristicsof the neural network and brain functional network features.

If the cerebral condition is determined to be abnormal, then thelocation of abnormal features (e.g., local or focal, lateralized, ordiffuse bilateral) and/or the severity of the abnormality (e.g., mild,moderate, or severe) may be identified in box 248. For example, thecondition may be identified as abnormal diffuse bilateral, abnormal lefthemisphere, or abnormal right hemisphere with slowing, seizures, and/oramplitude disturbance. An indication of the classification results maythen be generated in box 251 for rendering on the system display. Forexample, a graphic display of the original EEG signal, local signalproperties, inter-hemispheric asymmetries, local network features,and/or global network features may be generated. A warning may begenerated when an abnormal condition has been indicated.

A summary (or report) of findings may be provided in several forms whichmay be selected by the user. For example, a default condition mayprovide a report labeled as normal or indicating the determined abnormalcategory classification (e.g., mildly abnormal, left hemisphere). Inaddition, a visual display of the anatomical location of theabnormalities may be provided graphically, using a color bar, grey scaleor other graphic display to indicate the severity of the abnormalities.Other graphical displays which provide maps of one or more individualsignal property may also be viewed. The results of the classificationsmay also be stored in box 251 for later access or retrieval to furtherevaluation, interpretations, and validation.

Referring to FIGS. 3A and 3B, shown are a flowchart illustratingexamples of feature extraction, network modeling, and classification ofthe flowchart of FIG. 2B. In FIG. 3A, EEG data is received for EEGsignal feature extraction in box 236. In box 303, local properties maybe determined from the EEG data including the univariate feature.Kaiser-Teager energy and power in the EEG may be computed for eachelectrode site and referenced to a common reference electrode (e.g.,average, Pz, other). Other implementations may use bipolar pairs andcalculate for each pair: Fp1-F3, F3-C3, C3-P3, P3-O1 and analogous forright side, Fp1-F7, F7-T3, T3-T5, T5, O1 and analogous for right side.

Local energy and power properties may be determined for each channel forcomparison to predetermined normative values for each channel.Abnormality of Teager energy and abnormalities for power would beexpected to be either higher or lower than the normative values.Abnormalities in the Teager energy to Power ration would be expected tobe lower than norms. Norms may be derived from EEG recordings obtainedfrom a normal test group, with appropriate age matching, or may be basedupon baseline recordings obtained in the same subject (in which case, achange from baseline would be detected).

For each EEG channel, the follow local property values may be computedin box 303:

-   -   Kaiser-Teager energy (KTE). KTE may be calculated for each        electrode channel. This value can be obtained for the entire        recorded frequency range as well as for each of the standard EEG        frequency bands (delta: 0-4 Hz, theta: 4-8 Hz, alpha: 8-13 Hz,        and beta: 13-30 Hz). The value for each channel may be compared        to normal values. If the values are outside of the normal range,        the degree of abnormality (e.g., based on standard deviations        (s.d.) from the mean) for each electrode channel can be        determined. The location (left cerebral hemisphere, right        cerebral hemisphere or bilateral) and degree of abnormality (1        s.d.≧x<2 s.d., 2 s.d.≦x<3 s.d, or x≧3 .s.d.) can be stored and        used in the final evaluation and report.    -   Power. Standard power for the entire frequency range (0 to 30        Hz) and for each of the standard EEG frequency bands (delta,        theta, alpha and beta) may be computed for each electrode        channel and compared to normal values for each respective        channel. If the values are outside of the normal range, the        degree of abnormality (based on standard deviations from the        mean) for each electrode channel can be determined. The location        (left cerebral hemisphere, right cerebral hemisphere or        bilateral) and degree of abnormality (1 s.d.≧x<2 s.d., 2        s.d.≦x<3 s.d, or x≧3 .s.d.) can be stored and used in the final        evaluation and report.    -   Pattern match regularity statistic (PMRS). One or more measures        of signal regularity or signal such as the PMRS may be generated        for the entire recorded frequency range as well as for each        standard EEG frequency band for each electrode channel and        compared to normal values. If the values are outside of the        normal range, the degree of abnormality (based on standard        deviations from the mean) for each electrode channel can be        determined. The location (left cerebral hemisphere, right        cerebral hemisphere or bilateral) and degree of abnormality (1        s.d.≧x<2 s.d., 2 s.d.≦x<3 s.d, or x≧3 .s.d.) can be stored and        used in the final evaluation and report.    -   Approximate entropy. One or more measures of signal entropy,        such as approximate entropy, may be measured for the entire        recorded EEG frequency range as well as for each of the standard        EEG frequency bands for each electrode channel and compared to        normal values. If the values are outside of the normal range,        the degree of abnormality (based on standard deviations from the        mean) for each electrode channel can be determined. The location        (left cerebral hemisphere, right cerebral hemisphere or        bilateral) and degree of abnormality (1 s.d.≧x<2 s.d., 2        s.d.≦x<3 s.d, or x≧3 .s.d.) can be stored and used in the final        evaluation and report.

Teager energy/power ratios are generated for each channel for entirefrequency range between 1 and 30 Hz in box 306. For each EEG channel,the KTE to power ratio may be calculated for the entire recordedfrequency range as well as for the standard EEG frequency bands iscalculated for each channel and compared to normal values. If the valuesare outside of the normal range, the degree of abnormality (based onstandard deviations from the mean) for each electrode channel candetermined. The location (left cerebral hemisphere, right cerebralhemisphere or bilateral) and degree of abnormality (1 s.d.≧x<2 s.d., 2s.d.≦x<3 s.d, or x≧3 .s.d.) can be stored and used in the finalevaluation and report.

Left-right univariate ratios may then be determined in box 309.Inter-hemispheric symmetry computation may be based upon univariatefeatures. Each of the quantitated measures of signal properties, such asthose described above, will be examined for inter-hemispheric symmetryby calculating the ration of the value obtained for each of theelectrode channels recorded from the left cerebral hemisphere to thesame value obtained for the homologous electrode channel in the righthemisphere.

In box 312, bivariate features between all channels may also bedetermined from the EEG data. Inter-hemispheric symmetry computationsmay be based upon bivariate features. Bivariate measures can be used toevaluate the relationship of signals obtained from each electrodechannel from the left cerebral hemisphere to that of the homologouselectrodes from the right cerebral hemisphere. These measures includemutual information, linear or nonlinear correlation, coherence, phaselocking index and phase lag index. The analysis can be made for theentire range of recorded frequencies as well as for each standard EEGfrequency band.

A network model may then be generated in box 315 based upon thebivariate feature values from box 312. A network model can be generatedas a weighted graph, based on one or more of the bivariate measuresrelating signal properties between each pair of electrodes, such that ameasure is generated for each electrode site (node) and all otherelectrodes in the recording. The weighted graph can be converted to abinary graph depicting the node pairs with the strongest association, asdefined by one or more bivariate measure, using a threshold in box 318.For example, a threshold of 0.75 may be used, such that the resultantbinary graph includes 25% of the total electrode pairs; in the case of afull set of electrodes, excluding midline electrodes, as defined by theInternational 10-20 System of electrode placement, the total number ofpairs is 171 and 43 pairs would be selected for the binary networkgraph.

In box 321, global network characteristics of the binary and/or weightednetwork graphs may be determined. These characteristics include, e.g.,clustering coefficient and minimum path length. These globalcharacteristic values are compared to norms in box 324 to determinewhether or not they are within the normal range. If the values areoutside of the normal range, the degree of abnormality (based onstandard deviations from the mean) for each electrode channel isdetermined. The location (left cerebral hemisphere, right cerebralhemisphere or bilateral) and degree of abnormality (1 s.d.≧x<2 s.d., 2s.d.≦x<3 s.d, or x≧3 .s.d.) will be stored and used in the finalevaluation and report.

In addition, characteristics for each node (or electrode) can bedefined, based on the following characteristics for each node: degree,path length to contralateral homologous electrode, and connectionstrength with contralateral homologous electrode. For each electrode,the degree of that node can be compared to the degree for the sameelectrode (node) in the normative dataset in box 327 and the location ofnodes whose properties do not match those of the normative dataset canbe determined.

In box 330, hubs of the binary and/or weighted network graphs may beidentified using one or more criteria for defining hubs, such as degree,betweeness, closeness, and eigen vector centrality. Electrodes whichexceed thresholds of the values for each respective measure can bedefined as a network hub. Network hubs identified in the recording canbe compared to a list of hubs obtained from a normative comparisondataset. Hubs present in the subject network which do not correspond tohubs in the normal datasets may be identified and their locationdetermined to be lateralized to one cerebral hemisphere, localizedwithin one cerebral, or present bilaterally. In addition, nodes in therecorded data which are not present in the normative datasets can beidentified and localized. In a similar fashion, the path length betweeneach electrode (node) and the homologous node in the contralateralhemisphere may be calculated. Values for each channel pair can becompared to those of the normal dataset and the location of those pairswhich differ significantly from the normal datasets can be determined.

The cerebral condition can be classified in box 245 of FIG. 3B basedupon comparison of values from boxes 309 and 312 of FIG. 3A. Aftersignals have been analyzed and results stored, the EEG data may beclassified, based on the composite results from all, or a subset, ofeach individual analysis: (1) local univiariate signal properties, (2)inter-hemispheric symmetry computations based on local univariate signalproperties, (3) inter-hemispheric symmetry based on computations ofbiviariate signal properties, local network properties, and globalnetwork properties. If all of these analyses are within acceptable rangeof normal values, the recording may be classified as normal.

In box 333, the bivariate feature of homologous pairs is compared withnormal values obtained from the same electrode pairs. If the values areoutside of the normal range, the degree of abnormality (based onstandard deviations from the mean) for each electrode channel may bedetermined. The location (left cerebral hemisphere, right cerebralhemisphere or bilateral) and degree of abnormality (1 s.d.≧x<2 s.d., 2s.d.≦x<3 s.d, or x≧3 .s.d.) can be stored and used in the finalevaluation and report. In box 336 of FIG. 3B, univariate ratios fromeach channel pair will be compared to normal ratios for that channelpair. If the ratios are outside of the normal range, the degree ofabnormality (based on standard deviations from the mean) for eachelectrode channel may be determined. The location (left cerebralhemisphere, right cerebral hemisphere or bilateral) and degree ofabnormality (1 s.d.≧x<2 s.d., 2 s.d.≦x<3 s.d, or x≧3 .s.d.) can bestored and used in the final evaluation and report.

If there are abnormalities identified, the EEG data will be classifiedas abnormal and assigned to one of several abnormal categories in box248 of FIG. 3B. For example, the abnormal categories may be defined asfollows: (1) left hemisphere abnormality, (2) right hemisphereabnormality, (3) bilateral abnormalities. EEG data with bilateralabnormalities may be further subclassified as (3a) bilateral symmetricalabnormalities, (3b) bilateral abnormalities, left greater than right and(3b) bilateral abnormalities right greater than left. In otherimplementations, the location of the abnormality may be made moreprecise, indicating the region(s) within the cerebral hemisphere(s)containing the abnormalities (e.g., left front-temporal abnormality).Abnormal EEG data may be further categorized as to the degree ofabnormality, based upon a weighted magnitude of property deviations fromnormal values as well as the number of properties which deviatesignificantly from normal values.

In box 342, regions of altered symmetry are identified and the severitydetermined based upon values from box 333. In addition, the severity ofthe abnormality is determined in box 339 based upon values for boxes324, 327, and 330 of FIG. 3A and box 336 of FIG. 3B. As discussed above,the degree of abnormality may be based on standard deviations from themean the values for each electrode channel. In some cases, a weighedcombination of the standard deviations may be used to indicate thedegree of abnormality. The results of the abnormality identification maybe stored in a data store or memory and may be used to generate anindication for the operator of the system 100.

The system 100 may be used for, but is not limited to, neurologicalassessment of common neurological presentations such as, e.g., acuteencephalopathies, subacute encephalopathies, focal lesions, ischemicevents, and chronic encephalopathies. Acute and subacuteencephalopathies include such disorders as those due to traumatic braininjuries, toxic encephalopathies (e.g. drug or alcohol toxicity),metabolic disorders (e.g. hypoglycemia, hyperglycemia, ketoacidosis,renal failure, hepatic failure, hypoxia, hypercapnea), acute or subacuteinfections of the brain (such as meningitis, encephalitis, and brainabscess), seizures, status epilepticus, stroke, transient ischemicattacks, and autoimmune disorders affecting the central nervous system.It may also be used for detecting mild disturbances of brain functionincluding, e.g., concussion following after head injuries. Informationobtained through the system 100 may be used to refine the differentialdiagnosis, formulate further workup (e.g. imaging procedures) andtreatment, and for purposes of triage and referral to appropriatefacilities and specialists. Other potential uses include brainmonitoring to evaluate level of alertness, screening of chronic cerebraldisorders such as, e.g., Alzheimer's disease and other chronicdementias, and assessment of excess daytime sleepiness and sleepdisorders.

One embodiment, among others, can be an altered mental status evaluator(AMSE). A unique feature of this application is its utility as a tool toassist physicians in the differential diagnosis of subjects in thehospital with acute unexplained persistent altered mental status (AMS).These subjects are commonly seen in the Emergency Room (ER), IntensiveCare Unit (ICU) or hospital wards. AMSE provides reliable identificationof EEG abnormalities that cause altered mental status includingsubclinical seizures, diffuse EEG slowing reflecting diffuseencephalopathy and focal EEG slowing reflecting focal brain dysfunction.AMSE then generates a report to assist the physician in diagnosis ofsubjects with altered mental status. These results can rapidly point thephysician to the differential diagnostic areas that should receive themost consideration initially (but should not preclude other avenues ofinvestigation). The physician correlates the results with their clinicalexamination and results of other studies to reach a final accuratediagnosis more quickly.

With reference to FIG. 4, shown is a schematic block diagram of aprocessor system 400 in accordance with various embodiments of thepresent disclosure. The processor system 400 includes at least oneprocessor circuit, for example, having a processor 403 and a memory 406,both of which are coupled to a local interface 409. To this end, theprocessor system 400 may comprise, for example, at least one computer orlike device. The local interface 409 may comprise, for example, a databus with an accompanying address/control bus or other bus structure ascan be appreciated. In addition, the processor system 400 includesoperator interface devices such as, e.g., a display device 412, akeyboard 415, and/or a mouse 418. In some implementations, the operatorinterface device may be interactive display 421 (e.g., a touch screen)that provides various functionality for operator interaction with theprocessor system 400. Various sensors such as, e.g., EEG electrodes 424may also interface with the processor system 400 to allow foracquisition of EEG signals from a subject. In some embodiments, the EEGelectrodes 424 may be an array of electrodes configured to be positionedabout the subject's head.

Stored in the memory 406 are both data and several components that areexecutable by the processor 403. In particular, stored in the memory 406and executable by the processor 403 are various application modules 427such as, e.g., an electrode application module 103, an EEG recordingmodule 106, a signal conditioning module 109, a signal analysis module112, and a condition classification module 115 of FIG. 1, and/or otherapplications. Also stored in the memory 406 may be a data store 430 andother data. In addition, an operating system 433 may be stored in thememory 406 and executable by the processor 403.

It is understood that there may be other applications that are stored inthe memory 406 and are executable by the processor 403 as can beappreciated. Where any component discussed herein is implemented in theform of software, any one of a number of programming languages may beemployed such as, for example, C, C++, C#, Objective C, Java®,JavaScript®, Perl, PHP, Visual Basic®, Python®, Ruby, Delphi®, Flash®,or other programming languages.

A number of software components are stored in the memory 406 and areexecutable by the processor 403. In this respect, the term “executable”means a program file that is in a form that can ultimately be run by theprocessor 403. Examples of executable programs may be, for example, acompiled program that can be translated into machine code in a formatthat can be loaded into a random access portion of the memory 406 andrun by the processor 403, source code that may be expressed in properformat such as object code that is capable of being loaded into a randomaccess portion of the memory 406 and executed by the processor 403, orsource code that may be interpreted by another executable program togenerate instructions in a random access portion of the memory 406 to beexecuted by the processor 403, etc. An executable program may be storedin any portion or component of the memory 406 including, for example,random access memory (RAM), read-only memory (ROM), hard drive,solid-state drive, USB flash drive, memory card, optical disc such ascompact disc (CD) or digital versatile disc (DVD), floppy disk, magnetictape, or other memory components.

The memory 406 is defined herein as including both volatile andnonvolatile memory and data storage components. Volatile components arethose that do not retain data values upon loss of power. Nonvolatilecomponents are those that retain data upon a loss of power. Thus, thememory 406 may comprise, for example, random access memory (RAM),read-only memory (ROM), hard disk drives, solid-state drives, USB flashdrives, memory cards accessed via a memory card reader, floppy disksaccessed via an associated floppy disk drive, optical discs accessed viaan optical disc drive, magnetic tapes accessed via an appropriate tapedrive, and/or other memory components, or a combination of any two ormore of these memory components. In addition, the RAM may comprise, forexample, static random access memory (SRAM), dynamic random accessmemory (DRAM), or magnetic random access memory (MRAM) and other suchdevices. The ROM may comprise, for example, a programmable read-onlymemory (PROM), an erasable programmable read-only memory (EPROM), anelectrically erasable programmable read-only memory (EEPROM), or otherlike memory device.

Also, the processor 403 may represent multiple processors 403 and thememory 406 may represent multiple memories 406 that operate in parallelprocessing circuits, respectively. In such a case, the local interface409 may be an appropriate network that facilitates communication betweenany two of the multiple processors 403, between any processor 403 andany of the memories 406, or between any two of the memories 406, etc.The local interface 409 may comprise additional systems designed tocoordinate this communication, including, for example, performing loadbalancing. The processor 403 may be of electrical or of some otheravailable construction.

Although the electrode application module 103, the EEG recording module106, the signal conditioning module 109, the signal analysis module 112,the condition classification module 115, and other various systemsdescribed herein may be embodied in software or code executed by generalpurpose hardware as discussed above, as an alternative the same may alsobe embodied in dedicated hardware or a combination of software/generalpurpose hardware and dedicated hardware. If embodied in dedicatedhardware, each can be implemented as a circuit or state machine thatemploys any one of or a combination of a number of technologies. Thesetechnologies may include, but are not limited to, discrete logiccircuits having logic gates for implementing various logic functionsupon an application of one or more data signals, application specificintegrated circuits having appropriate logic gates, or other components,etc. Such technologies are generally well known by those skilled in theart and, consequently, are not described in detail herein.

Although the flowcharts of FIGS. 2A-2B and 3A-3B show a specific orderof execution, it is understood that the order of execution may differfrom that which is depicted. For example, the order of execution of twoor more blocks may be scrambled relative to the order shown. Also, twoor more blocks shown in succession in FIGS. 2A-2B and 3A-3B may beexecuted concurrently or with partial concurrence. Further, in someembodiments, one or more of the blocks shown in FIGS. 2A-2B and 3A-3Bmay be skipped or omitted. In addition, any number of counters, statevariables, warning semaphores, or messages might be added to the logicalflow described herein, for purposes of enhanced utility, accounting,performance measurement, or providing troubleshooting aids, etc. It isunderstood that all such variations are within the scope of the presentdisclosure.

Also, any logic or application described herein, including the electrodeapplication module 103, the EEG recording module 106, the signalconditioning module 109, the signal analysis module 112, the conditionclassification module 115, and/or application(s), that comprisessoftware or code can be embodied in any non-transitory computer-readablemedium for use by or in connection with an instruction execution systemsuch as, for example, a processor 403 in a computer system or othersystem. In this sense, the logic may comprise, for example, statementsincluding instructions and declarations that can be fetched from thecomputer-readable medium and executed by the instruction executionsystem. In the context of the present disclosure, a “computer-readablemedium” can be any medium that can contain, store, or maintain the logicor application described herein for use by or in connection with theinstruction execution system. The computer-readable medium can compriseany one of many physical media such as, for example, magnetic, optical,or semiconductor media. More specific examples of a suitablecomputer-readable medium would include, but are not limited to, magnetictapes, magnetic floppy diskettes, magnetic hard drives, memory cards,solid-state drives, USB flash drives, or optical discs. Also, thecomputer-readable medium may be a random access memory (RAM) including,for example, static random access memory (SRAM) and dynamic randomaccess memory (DRAM), or magnetic random access memory (MRAM). Inaddition, the computer-readable medium may be a read-only memory (ROM),a programmable read-only memory (PROM), an erasable programmableread-only memory (EPROM), an electrically erasable programmableread-only memory (EEPROM), or other type of memory device.

It should be emphasized that the above-described embodiments of thepresent disclosure are merely possible examples of implementations setforth for a clear understanding of the principles of the disclosure.Many variations and modifications may be made to the above-describedembodiment(s) without departing substantially from the spirit andprinciples of the disclosure. All such modifications and variations areintended to be included herein within the scope of this disclosure andprotected by the following claims.

EEG Analysis of Brain Dynamical Behavior with Applications in EpilepsyLIST OF ABBREVIATIONS

-   AED Anti-epileptic drug-   ANOVA Analysis of variance-   AR Autoregressive-   ARMA Autoregressive moving average-   ARIMA Autoregressive integrated moving average-   CC Cross-correlation-   CPS Complex partial seizure-   EEG Electroencephalogram-   EMU Epilepsy monitoring unit-   EWS Evolutionary wavelet spectrum-   FN False negative-   FPR False positive rate-   LSW Local stationary wavelet-   GABA Gamma-Aminobutyric acid-   MSC Mean square coherence-   MTLE Mesial temporal lobe epilepsy-   PLAI Phase lag index-   PMRS Pattern-match regularity statistics-   PNES Psychogenic nonepileptic seizures-   STLmax Short-term maximum Lyapunov exponent-   SVM Support vector machine-   TN True negative-   TP True positive

Chapter 1 Introduction Electroencephalogram

Electroencephalogram (EEG) is a technology measuring the voltage ofactivation from a set of neurons within the brain. EEG electrodes canrecord brainwaves from either the scalp or, through invasive procedures,deeper layers of the brain tissue. Depending on the location and thetype of the electrodes, EEG signals can reveal different levels ofneuronal activities. Scalp EEG is the least invasive methodology and is,therefore, widely used in many fields of neuroscience includingpsychology, epilepsy, brain machine interface, and sleep research. EEGhas become a more popular means of recording interesting brain activityover the last few decades because improvements in computing power andstorage capacity make possible sophisticated analyses of very largeamounts of EEG data.

Like electrocardiogram (ECG), another widely used bioelectrical signalin monitoring and diagnosing heart abnormalities, numerous electrodescan be applied on one subject so that the EEG data have both temporaland spatial information of brain activity. The location of electrodeplacement on a scalp usually follows the 10-20 system. Distinctlandmarks of a head are first identified and then electrodes are placedat 10% or 20% distance intervals along the landmarks. Through analyzingthe multi-channel EEG data of a short period of time, a temporary brainnetwork can be identified. Analyzing long-term EEG recording data canreveal the dynamics of the brain network as the subject goes throughmany states of mental activity or pathological manifestation.

The EEG data included in this disclosure contain both scalp andintracranial EEGs. All of them are continuous long-term recordings fromeither Allegheny General Hospital (AGH, Pittsburgh, Pa.) or the MedicalUniversity of South Carolina (MUSC, Charleston, S.C.). Depending on thestudy design, long-term continuous or extracted EEG recordings wereanalyzed. All scalp EEG data were recorded using 19 channels andintracranial EEG data used respectively different channels covering thedifferent brain areas depending on the specific case of each epilepticpatient. All recordings are from human beings who are older than 18years old and experienced epileptic or psychological non-epilepticseizures.

Epilepsy

Epilepsy is a brain disorder, instead of a disease, characterized mostlyby recurrent and unforeseen interruptions of normal brain function. Ithas been known from ancient times and was, at one point, attributed todivine intervention until the great Greek physician Hippocrates realizedthat it was a disorder of the brain. It is also the second most commondisorder of the central nervous system and affects about 0.4-1% of thepopulation.

The sudden interruption of brain function is termed a seizure. Manyother physical or psychological sudden and severe events are also calledseizures, such as a heart seizure. The meaning of a seizure, going backto its' origin in Greek, is “to take hold”. An epileptic seizure is themanifestation of epilepsy and is due to abnormally excessive orsynchronous neuronal activity in the brain. The duration of themanifestation of an epileptic seizure is sometimes called an ictalperiod. Most epileptic seizures can be recorded by the electrodes of anEEG and shown in the EEG data. Through reading and analyzing theepileptic EEG data, more information is gained about the onset patternand type of epilepsy. During the onset of an epileptic seizure (ictalperiod), one may observe that the synchronous and rhythmic dischargesoriginate from one part of the brain (partial, focal orlocalization-related seizures) or begin simultaneously in both sides ofthe hemispheres (generalized seizures). Focal seizures after the onsetmay remain localized within one part of the brain or propagate to theother side of the hemisphere and cause a wider range of synchronousneuronal activity (secondarily generalized seizures). If a seizure isconfined to a limited area of the brain, it usually causes relativelymild and transient cognitive, psychic, sensory, motor or autonomicimpairment. However, generalized seizures cause altered consciousnessand accompany a variety of motor symptoms from the jerking of limbs tostiff or convulsive movements of the whole body. During the ictal periodof focal seizures, consciousness may be unaffected (simple partialseizures) or be impaired (complex partial seizures). Few patients withfocal seizures feel unusual sensations (auras) when experiencing thebeginning stage of an ictal period. Most seizures occur in a nearlyunpredictable manner, stressing the patients as well as the peoplearound them.

Since an epileptic seizure is characterized by excessively synchronousneuronal activity, researchers hypothesize that the dysfunction of theinhibition mechanism of neurons favors the development of seizures.Gamma-aminobutyric acid (GABA) is the most prevalent inhibitoryneurotransmitter in human brains. GABAergic inhibition is one of themain control mechanisms influencing the excitability of neurons. Manycompounds that increase the GABAergic inhibition are used asanti-epileptic drugs (AED). Although many studies confirm therelationship between GABAergic inhibition and the development ofepilepsy, about 20-30% of epileptic patients failed to benefit from theadministration of AED. These latter patients expressed poor quality oflife due to both depression and neurotoxicity from AED. Althoughpatients have different seizure onset frequencies and severities,seizure types and frequencies did not correlate with the score ofquality of life. These study results suggest that the frequency orseverity of seizures is not the primary cause of patient depression, butrather the sense of uncertainty regarding seizure onsets and theconstant fear of unpredictable onset. The presence of correct seizureonset warnings should alleviate the degree of depression and increasequality of life by decreasing the amount of uncertainty and offeringpatients the ability to prepare for seizure onset.

Chapter 2 The Properties and Analysis of Non-Stationary Time SeriesForeword

Most phenomena happening with time can be measured and recorded as timeseries, such as wind speed, heart rate, stock exchange price, and EEGsignals. By analyzing a time series, insight may be gained into an arrayof interesting and complex phenomena. A human brain is a complex systemthat responds simultaneously to many external inputs as well as numerousinternal processes. Analysis and quantification of the characteristicsof EEG signals offer a window to observe the complex interactionsamongst neurons in a brain. For example, in awake and asleep EEGsignals, one can recognize that the signals usually have swiftlychanging backgrounds. Therefore, EEG signals are usually regarded asnon-stationary time series consisting of the many distinct backgroundactivities that occur when a subject experiences different phases ofpsychological events or physiological states. Analyzing and forecastingnon-stationary time series have been challenging for researchers. Manypossible methodologies have been developed to alleviate the difficultiesof analyzing a non-stationary process. Chapter 2 introduces the basicproperties and analysis of a non-stationary times series. Since therapid growth in the power of computational machines, the ability toinvestigate very data rich EEG signals with sophisticated andcomputation-demanding techniques has been gained.

For an interesting time series, it is convenient to assume that thereare exact mathematical models generating the time series. Understandingthe underlying mechanism of these processes is very crucial inengineering and science. It is also intuitive to presume that there arerules governing the route of a time series and once those rules areunderstood, a forecast of how the time series may develop may beimproved. However, real life phenomena are rarely that simple andusually involve many parameters, noise, and unobserved variables thatcan influence the time series with different time lags to a differentextent. Still, with careful control of the recording environment, it isfeasible to reduce many unobserved covariates and try to model the timeseries as a stochastic process. The more precise the time series model,the smaller the forecasting error that can be achieved.

To understand all the principles and rules controlling these black-boxsystems at once is nearly impossible. Therefore, simple non-stationaryparadigms are introduced before moving to the more complicated timeseries. As aforementioned, a straightforward assumption is that aprocess possesses certain fixed properties in itself so that it evolveswith principles. One of the most convenient properties in a stochasticprocess is stationarity. Stationarity in the strict sense means thosestatistics of a stochastic process do not change with time as describedin Equation 2-1.

ƒ_(x) ₁ _(,x) ₂ _(, . . . ,x) _(n) (t ₁ ,t ₂ , . . . ,t _(n))=ƒ_(x) ₁_(,x) ₂ _(, . . . ,x) _(n) (t ₁ +τ,t ₂ +τ, . . . ,t _(n)+τ)  (2-1)

In Equation 2-1, ƒ denotes the probability mass function and T denotes atime lag of any size. When a process fails to fulfill Equation 2-1, theprocess is considered not stationary in the strict sense.

Stationary processes possess relatively stable statistical propertiescompared to non-stationary ones; such as a fixed mean and variance.Nevertheless fixed mean and variance by themselves do not implystationarity in the strict sense. With those fixed statistics, it ispossible to begin predicting the foregoing process with more confidence.The first section of Chapter 2 focuses on some background for methodsthat will be introduced later. In the second section of Chapter 2, twowell-known practical methods to forecast non-stationary processes arepresented.

Properties of Time Series Stationarity

Stationarity in the strict sense is almost impossible in real life. Evena parametric simulated process can have stationarity only theoretically,not to mention a time series recorded from real world signals which havemuch more complex mechanisms. Although one can hardly have a stationaryprocess in the strict sense, the concept of stationarity may be utilizedwhen studying non-stationary time series because of its convenientstatistical properties. A process may be claimed to have stationarity ina wide sense if it has not only a fixed mean and variance but alsotime-invariant auto-covariance. One reason for keeping the concept ofstationarity in a stochastic process is that a parametric model canhardly include all variables that affect the process, not to mentionsome factors that are beyond the researcher's control. To take alleffects into account, a stochastic process should be modeled withseveral deterministic and random terms. One way to implement this is touse the concept of piecewise stationarity. Even though the whole processis non-stationary, each segment of a reasonable length may be treated asif they were stationary, after disintegrating it into segments. Based ona piecewise stationary concept in the time domain, one can make ananalogue to the time and frequency domain, which is the idea of localstationary wavelet analysis which will be introduced in thenon-parametric method section.

Local Stationarity

One of the crucial parameters for using the concept of piecewisestationarity is how to choose the length of a segment that is neithertoo short for calculating statistical estimates nor too long so that theprocess is no longer stationary. It should be noted that the criterionfor selecting the segment length should always be dependent on themethod used in order to model the non-stationary process. For example,for a quasi-stationary process, the underlying assumption is that theparametric properties might change their state between different timesegments. As a result, an ideal segment should not be too long tooverlap segments having different parametric values. Compared toquasi-stationary processes, parameters of a time varying autoregressivemodel do not suddenly change, but gradually evolve along with time.Proper selection of a method to model a process depends on some basicproperties of the raw data. For instance, a seismic wave recording or anelectroencephalograph recording containing a seizure shows a suddenstate change. As a result, it is reasonable to use a quasi-stationarymethod for modeling. For other processes, such as the population of somespecies or the price of a stock, one might only be interested in theirsmooth evolution property instead of their sudden change due to specialor dramatic events, so it would be reasonable to use a time-varyingautoregressive model that will be introduced later. Those methods ofmodeling a non-stationary process are called parametric methods. Inaddition, one can use non-parametric methods which do not assume theexistence of a model or distribution family in a process. Compared toparametric methods, non-parametric methods still involve parameterchoosing to some extent but not as much as parametric methods. If oneintends to observe a non-stationary process as a piecewise, quasi, orevolutionary stationary process, a very well-known and widely usednon-parametric tool called time-dependent spectrum in the time-frequencyanalysis field can be used.

Parametric Methods for Analyzing Non-Stationary Time Series

The benefits of using a parametric method are its simplicity and lucidstructure. The main feature of parametric methods is the use ofparameters as structural elements to introduce interpretation andinteraction so that one can reduce ambiguity within the data itself.This, however, might be achieved at the risk of misinterpretation. Someother parametric approaches are based on specific assumptions thatprocess outputs belong to the family or families of the distributionsused in the model or can be described as the output of a stochasticprocess, whereas non-parametric approaches do not impose specificconstraints with regards to the distribution family. One has to firstcarefully scrutinize the statistical properties of the observed data andthe assumption for the parametric method being applied for estimation.In sum, parametric methods are beneficial only if there is a strongreason to look at the process in a particular structural way or as amember of a certain distribution family. As one can imagine, aparametric approach might not perform well on a process having a hugefuture random term that does not constantly fall into any distributionfamily and has a weak link between current and future states of theprocess. Time-Varying Autoregressive (AR) models, just like a stationaryprocess, can use parametric or non-parametric means to start theanalysis on a non-stationary process. One of the most commonly usedmodels for simulating or fitting non-stationary time series is thetime-varying AR model. A p^(th) order linear time-varying AR model canbe written in a general form as Equation 2-2.

$\begin{matrix}{x_{k} = {{\sum\limits_{i = 1}^{p}\; {\varphi_{i,k}x_{k - i}}} + ɛ_{k}}} & \left( {2\text{-}2} \right)\end{matrix}$

In Equation 2-2, φ_(i,k) are time-varying parameters that differentiatea time-varying AR model from an AR model that has a fixed coefficientφ_(i) instead of φ_(i,k). The main difficulty in using a parametricmodel for analysis lies in choosing proper values for the parameters. AsEquation 2-2 shows, p and φ_(i,k) are parameters which need to be chosenbased on the training data while ε_(k) is a zero-mean white noiseprocess. One advantage of a parametric method is that it renders a luciddynamic relationship between data points and gives some sense of howdata evolve over time. However, the model may not be unique because ofthe uncertainty of parameters. One should bear in mind that models ofdifferent orders may be fit into the same datasets. In this case thetime varying AR model parameter values (φ_(i,k)) might be significantlydifferent especially for a complex dataset. Several methods have beenproposed to find the optimal order of an AR model such as Akaikeinformation criterion, Bayesian information criterion, and thelikelihood ratio test. Once the order of the AR model is decided, φ_(i)are calculated from Yule-Walker equations. Other methods exist forestimating regression coefficients, such as Burg's method, least-squaresapproach, modified covariance method, covariance method, parametricspectral estimation method, Prony's method, etc. However, choosing theorder of a time-varying AR model is more complex than that of an ARmodel. A recently proposed method for time-varying AR model orderuncertainty is the generalized likelihood ratio test.

Autoregressive Integrated Moving Average (ARIMA) Model:

Modeling an interesting time series as a non-stationary processes whosederivative, of proper order, is a stationary process. An autoregressivemoving average model (ARMA) can be represented in Equation 2-3.

$\begin{matrix}{{x_{n} = {{- {\sum\limits_{k = 1}^{p}\; {a_{k}x_{n - k}}}} + {\sum\limits_{k = 0}^{q}\; {b_{k}u_{n - k}}}}},\mspace{14mu} p,{q \in Z}} & \left( {2\text{-}3} \right)\end{matrix}$

In Equation 2-3, {u_(n)} is a Gaussian white noise process. An ARIMA canexhibit not only ARMA but also non-stationarity having neither fixedmean nor fixed variance. A very important property of ARIMA is that ifone takes the d^(th) derivative of the process it will result in astationary process that can be represented by another ARMA model{x_(n)}. This relationship can be represented in Equation 2-4.

$\begin{matrix}{{{\Delta^{d}z_{n}}\overset{\Delta}{=}{x_{n} = {{- {\sum\limits_{k = 1}^{p}\; {a_{k}x_{n - k}}}} + {\sum\limits_{k = 0}^{q}\; {b_{k}u_{n - k}}}}}},\mspace{14mu} p,q,{d \in Z}} & \left( {2\text{-}4} \right)\end{matrix}$

This relationship renders a powerful linkage for changing anon-stationary process into a stationary process which can be moreeasily analyzed. However, one should notice that the high frequencycomponents of the process would be amplified after differentiating.Considering an ARMA (p,q), with q≦p, it is possible that the AR part hasmore terms than the MA part. The realization of an ARMA state spacemodel is given in Equation 2-5.

$\begin{matrix}{{{\overset{\sim}{x}}_{k + 1} = {{F_{k}{\overset{\sim}{x}}_{k}} + {\begin{bmatrix}{b_{1,k} - a_{1,k}} \\{b_{2,k} - a_{2,k}} \\\vdots \\{b_{p,k} - a_{p,k}}\end{bmatrix}{\overset{\sim}{u}}_{k}}}},} & \left( {2\text{-}5} \right) \\{{z_{k} = {{\begin{bmatrix}1 & 0 & \ldots & 0\end{bmatrix}{\overset{\sim}{x}}_{k}} + u_{k}}},} & \left( {2\text{-}6} \right)\end{matrix}$

In Equation 2-5, {tilde over (x)}_(k)={x_(k), x_(k−1), . . . ,x_(k−p+1)}^(T) and {ũ_(k)}, {u_(k)} are both Gaussian white noiseprocesses. Tilda indicates a vector. If the transition matrix is denotedas F_(k), the noise coupling matrix as G_(k) and the observation matrixas {tilde over (h)}=[1, 0, . . . 0]. Using Equation 2-5 and Equation2-6, Equation 2-3 and Equation 2-4 become Equation 2-7 and Equation 2-8.

{tilde over (x)} _(k+1) =F _(k) {tilde over (x)} _(k) +G _(k) ũ_(k),  (2-7)

z _(k) ={tilde over (h)}{tilde over (x)} _(k) +u _(k)  (2-8)

We can make a comparison between ARIMA and time-varying AR which has thedynamic matrix form as in Equations 2-9 and 2-10.

{tilde over (x)} _(k+1) =F _(k) {tilde over (x)} _(k)+ε_(k),  (2-9)

z _(k) ={tilde over (h)}{tilde over (x)} _(k),  (2-10)

In Equation 2-9, ε_(k) is a zero-mean innovation with time-varyingvariance and F_(k) is described in Equation 2-11.

$\begin{matrix}{F_{k} = {\begin{bmatrix}\varphi_{1,k} & \varphi_{2,k} & \ldots & \varphi_{{p - 1},k} & \varphi_{p,k} \\1 & 0 & \ldots & 0 & 0 \\0 & 1 & \ldots & 0 & 0 \\\vdots & 0 & \ddots & \vdots & \vdots \\0 & 0 & \ldots & 1 & 0\end{bmatrix}.}} & \left( {2\text{-}11} \right)\end{matrix}$

Non-Parametric Methods for Analyzing Non-Stationary Time Series

Non-parametric methods assume no specific distribution or fixedstructure existing for a process. Some parameter settings still exist inthe non-parametric methods but the parameters are mainly for thesensitivity aspects of the analysis rather than relating to wholeprocess. As a result, non-parametric methods are less affected by thechoice of assumptions or parameter value compared to a parametricmethod. In sum, non-parametric methods do not explain the mechanismbehind a time series and are more data-driven.

Wavelet Transforms

Wavelet analysis is one of the most utilized tools for time-frequencyanalysis. Most of the stationary processes are generated and recorded inthe time domain. However, a weak stationary process can be characterizedby its' frequency components. That is why Fourier transform is appliedin almost every field of science and engineering. The Fourier transformrepresents a process formulating a superposition of frequencies.

X(ƒ)=∫_(−∞) ^(∞) x(t)e ^(−2πift) dt  (2-12)

x(t)=∫_(−∞) ^(∞) X(ƒ)e ^(2πift) df  (2-13)

In Equation 2-12, i denotes an imaginary unit. Although therepresentation through Fourier transform can clearly show those weightedcomponents of a process in terms of frequencies, it does not include anytime component. One can hardly perceive when those frequency componentsparticipate in a process. Especially for real data such as seismic wavesor biomedical signals, a functional transform showing results in boththe time and frequency domain can be much more informative andintuitive. A very widely used time-frequency analysis is the wavelettransform which renders a perspective with both time and frequencyresolutions through imposing a mask function before the original signal.

X(t,ƒ)=∫_(−∞) ^(∞) w(t−τ)x(τ)e ^(−t2πft) dτ,  (2-14)

In Equation 2-14, w(t−τ) is a window function. When w(t−τ) is arectangular window the above equation is also called short-time Fouriertransform, windowed Fourier transform, or time-dependent Fouriertransform. One may observe the characteristic of the wavelet fromEquation 2-14 that it has a t parameter that can characterize the signalat a certain range on the time domain. A real life signal such asseismic waves show characteristics of a short impulsive feature that canhardly be described through Fourier transform decomposing the seismicwave into a combination of sinusoidal functions ranging from negativeinfinity to infinity. Because of the dual property of a wavelet in boththe time and frequency domain, one cannot use a wavelet fortime-frequency analysis only, but also for de-nosing a process that hasnoise with a certain unique frequency band and that occurs sporadicallyalong the time frame.

Local Stationary Wavelets

Assuming that a process {X,}, tεZ possesses stationarity, second-orderstatistics such as variance or auto-covariance may be estimated throughonly a segment of observation and have confidence in the estimatedstatistics because its stationarity has rendered characteristicquantitative homogeneity. Having longer term observation, statistics maybe estimated better through elongating the observational period. Incontrast, for a non-stationary process, assuming that the variance of aprocess changes slowly over time, saying that variance is a continuousand differentiable function of time, the variance could still beestimated with certain precision by gathering information around thetime point of interested variance. This demand exactly fits the propertyof wavelet transform which preserves the local characteristics of aprocess to some extent on both the time and frequency domainrespectively. A statistically slowly changing process can be called alocal stationary process. For a covariance-stationary process, X_(t),there is a Cramér representation as Equation 2-15.

$\begin{matrix}{{X_{t} = {\int_{({{- \pi},\pi})}{{A(\omega)}{\exp \left( {{\omega}\; t} \right)}\ {{Z(\omega)}}}}},{t \in Z}} & \left( {2\text{-}15} \right)\end{matrix}$

In Equation 2-15, exp means exponential and Z(ω) is a stochastic processwith orthonormal increments. Non-stationary processes can be writtenalmost the same as Equation 2-15 but A(ω) is replaced by a function oftime. Following this rationale, a non-stationary process can berepresented using local stationary wavelet and non-decimated discretewavelets Ψ_(j,k)(t), j=−1, −2, . . . , kεZ where j indicates the scaleand k indicates time location. To express a local stationary process,one can write it as Equation 2-16.

$\begin{matrix}{{X_{t} = {\sum\limits_{j = {- J}}^{- 1}\; {\sum\limits_{k \in Z}\; {w_{j,k}{\psi_{jk}(t)}\xi_{jk}}}}},} & \left( {2\text{-}16} \right)\end{matrix}$

In Equation 2-16, ξ_(jk) is the mutually orthogonal zero mean randominnovation. From Equation 2-15 and Equation 2-16, one should notice thatw_(j,k) in Equation 2-15 corresponds to A(ω) in Equation 2-15 and bothof them indicate the amplitude of each analysis component. One shouldalso notice that Ψ_(j,k)(t) is a replacement of the Fourier harmonicstem, exp(iωt), in Equation 2-15. The subscript indices j and k inw_(j,k) represent scale and time location respectively just like thosein Ψ_(j,k) (t). The reason for using a non-decimated discrete wavelet isthat it can be shifted to any time point and not only confined in shiftsof 2^(−j) compared with a discrete wavelet. Please note thatnon-decimated wavelets do not have orthogonalities and are anovercomplete collections of shifted vectors. An example of a discretenon-decimated wavelet is the Haar wavelet. Nason also proposed anevolutionary wavelet spectrum (EWS) to quantify how the size of w_(j,k)changes over time as in Equation 2-17.

S _(j)(z)=S _(j)(k/T)≈w _(j,k) ².  (2-17)

Please note that in EWS, S_(j)(z), still has resolution on both thefrequency and time domain but with a different time scale which rescalesthe whole observation time into zε(0,1). Rescaled time, z, is calculatedby dividing k by the whole observation time length, T, and k 0, 1, . . ., T−1. The goal is to estimate the variance of a local stationaryprocess. To accomplish that, one way is to put together the localizedinformation around the time point of interest, with the concept of localstationary wavelet and EWS, let c(z,τ) represent localizedautocovariance and be defined as Equation 2-18.

$\begin{matrix}{{\lim\limits_{T\rightarrow\infty}{{cov}\left( {X_{{\lbrack{zT}\rbrack} - \tau},X_{{\lbrack{zT}\rbrack} + \tau}} \right)}} = {{c\left( {z,\tau} \right)} = {\sum\limits_{j = {- \infty}}^{- 1}\; {{S_{j}(z)}{\Psi_{j}(\tau)}}}}} & \left( {2\text{-}18} \right)\end{matrix}$

In Equation 2-18, Ψ_(j)(τ) is the autocorrelation function of Ψ_(j,k)(t)and [.] denotes the integer part of the real number.

$\begin{matrix}{{\Psi_{j}(\tau)} = {\sum\limits_{k = {- \infty}}^{\infty}\; {{\psi_{jk}(0)}{\psi_{jk}(\tau)}}}} & \left( {2\text{-}19} \right)\end{matrix}$

Once localized autocovariance is defined, localized variance at time zcan be defined as described in Equation 2-20.

$\begin{matrix}{{v(z)} = {{c\left( {z,0} \right)} = {\sum\limits_{j = {- \infty}}^{- 1}\; {S_{j}(z)}}}} & \left( {2\text{-}20} \right)\end{matrix}$

Further evidence supporting that a time-frequency domain autocovariancein terms of wavelets preserves the intuitive sense of an autocovariancein the time domain can be seen through the relationship between Allanvariances, which is defined as Equation 2-21, and Haar waveletvariances.

$\begin{matrix}{{{{\overset{\_}{X}}_{t}(\tau)} = {\frac{1}{\tau}{\sum\limits_{n = 0}^{\tau - 1}\; X_{t - n}}}},} & \left( {2\text{-}21} \right)\end{matrix}$

It is worth noting that Allan variance involves merely localized timedomain representation. However, a Haar wavelet involves both localizedtime and frequency domain representations. In this case, an Allanvariance can be expressed as in Equation 2-22.

$\begin{matrix}{{{\sigma_{X}^{2}(\tau)} = {\frac{1}{2}{E\left( {{{{\overset{\_}{X}}_{t}(\tau)} - {{\overset{\_}{X}}_{t - \tau}(\tau)}}}^{2} \right)}}},} & \left( {2\text{-}22} \right)\end{matrix}$

A Haar coefficient on the finest scale or smallest dilation, j=−1, attranslation or location, can be expressed as Equation 2-23.

d _(−l,k)=1/[√{square root over (2)}(X _(2k+1) −X _(2k))], k=0, . . .,T/2−1  (2-23)

Then it can see that they have the relationship as Equation 2-24.

var{{circumflex over (d)} _(j,k) }=E{circumflex over (d)} _(j,k)²=τ_(j)σ² _(X)(τ_(j)),  (2-24)

Forecasting Non-Stationary Time Series

Based on the background mentioned above, prediction algorithms fornon-stationary processes may be presented. A Kalman filter can be usedfor prediction on a more parametric basis. Compared to a Kalman filter,local stationary wavelets do not necessarily need to fit data with aparametric model and may concern only measurement per se. In contrast toa parametric analysis, a non-parametric analysis forsakes the idea thatthere is a model governing the evolution of data and that, in any case,the observation should be allowed to characterize the data by itself.

Kalman Filter

The Kalman filter is a recursive filter that can estimate linear dynamicparameters given some noisy measurements. Here the word, filter, shouldnot be taken as only a passive data processing algorithm, but think ofit as an active computer program that gathers information from inputsand then optimally estimates the transition of a system. In the previoussubsection introducing ARIMA, a non-stationary process can be treated asa stationary ARMA after taking appropriate orders of derivative on theARIMA. As a result, first analyze only ARMA without loss of generalityand then transform the result back to the ARIMA model. If an ARIMA failsto be differentiated properly into an ARMA, one can still implement aKalman filter onto an ARIMA process even if the ARIMA has missingobservations. Since the Gaussianity and linearity possessed in a Kalmanfilter are independent from the non-stationarity of a process, linearityand Gaussianity in the Kalman filter can be assumed and used to analyzea non-stationarity process.

According to work by McGonigal and Ionescu in 1995, a state space modelshould be found first before applying a Kalman filter for prediction. Itwas proposed that there was a simplified model based on Masanao Aoki'smethod that can model an ARMA state space. In Huang et al, 1995, theauthors implemented a Kalman filter to predict a time-varying AR (2)process which is a non-stationary and Gaussian-Markovian process and canbe described in a general form as Equation 2-25.

$\begin{matrix}{{x_{k} = {{\sum\limits_{i = 1}^{p}\; {\varphi_{i,k}x_{k - i}}} + ɛ_{k} + \sigma_{k}}},} & \left( {2\text{-}25} \right)\end{matrix}$

If one compares Equation 2-25 to Equation 2-2, one should notice thatthe time-varying model mentioned in this section has an extra term,σ_(k). Huang added σ_(k), a time-varying term, to better track themoving trend in the data of wind speed which is one of the applicationsin his work. To describe this time-varying AR (2) model in a matrixform, first define a vector space Ψ(φ₁ φ₂ σ)^(T) containing parametersthat need to be decided later and then represent the differences ofparameters as they evolve with time using their dynamic relationship.

$\begin{matrix}{{\begin{pmatrix}\Psi_{k} \\\Gamma_{k}\end{pmatrix} = {{\begin{pmatrix}I & I \\0 & H\end{pmatrix}\begin{pmatrix}\Psi_{k - 1} \\\Gamma_{k - 1}\end{pmatrix}} + {\begin{pmatrix}0 \\I\end{pmatrix}\Omega_{k}}}},} & \left( {2\text{-}26} \right)\end{matrix}$

In Equation 2-26, I is the identity matrix, Ω is a zero-mean white noisevector. If H is defined as an identity matrix, Γ_(k) is a random walkprocess and then Ψ_(k) ends up being an integrated random walk process.If H is a zero matrix, then Γ_(k) should be a white Gaussian noiseprocess and Ψ_(k) ends up being a random walk process. Please note thatthe role of Γ_(k) is the first order difference of Ψ_(k) and itdescribes the changes between Ψ_(k) and Ψ_(k−1) as time moves on. Anon-stationary process can be predicted, which is a time-varying AR (2)model, once the crucial step of choosing parameters in the model is doneand this can be accomplished by Kalman algorithm. A general way ofimplementation of a Kalman filter starts with setting up the parametermodel in the state space form.

{tilde over (z)} _(k) =A{tilde over (z)} _(k−1) +B{tilde over (ζ)}_(k),  (2-27)

{tilde over (y)} _(k) =C{tilde over (z)} _(k)+{tilde over(η)}_(k),  (2-28)

In Equation 2-28,

$\begin{matrix}{{{\hat{X}}_{{t + h - 1},T} = {\sum\limits_{s = 0}^{t - 1}\; {b_{{t - 1 - s},T}^{({t,h})}X_{s,T}}}},} & \left( {2\text{-}33} \right)\end{matrix}$

is the state vector and y_(k) is the observation vector of dimension p,meaning {tilde over (y)}_(k)=(x_(k) . . . x_(k−p+1))^(T). Given theestimation of z_(k) as {circumflex over (z)}_(k), one can write theestimation dynamics in terms of a current estimation, a past estimation,and a dynamic transform matrix A.

{circumflex over (z)} _(k|k−1) =A{tilde over (z)} _(k−1),  (2-29)

Then define the covariance matrix of estimation errors P_(k)=

({circumflex over (z)}_(k)−{tilde over (z)}_(k))^(T)({circumflex over(z)}_(k)−{tilde over (z)}_(k))

, state disturbance matrix R=

{tilde over (ζ)}_(k) ^(T){tilde over (ζ)}_(k)

, where

denotes expectation, and the noise variance matrix of measurement Q=

{tilde over (η)}_(k) ^(T){tilde over (η)}_(k)

. The updated state vector estimation, {circumflex over (z)}_(k), can becalculated through Young's recursive algorithm

{circumflex over (z)} _(k) ={circumflex over (z)} _(k−1) +P _(k|k−1) C^(T) [CP _(k|k−1) C ^(T) +R] ⁻¹ [y _(k) −C{circumflex over (z)}_(k|k−1)],  (2-30)

P _(k) =P _(k|k−1) −P _(k|k−1) C ^(T) [CP _(k|k−1) C ^(T) +R] ⁻¹ CP_(k|k−1),  (2-31)

P _(k|k−1) =AP _(k|k−1) A ^(T) +BQB ^(T),  (2-32)

Local Stationary Wavelet (LSW)

Given a series of t measurements, in which it is desired to use those tmeasurements to predict the process h step after t, an interestedprocess can be written as X_(0,T), . . . , X_(t−1,T) where T=t+h andthen define a linear predictor using Equation 2-33.

$\begin{matrix}{{{\hat{X}}_{{t + h - 1},T} = {\sum\limits_{s = 0}^{t - 1}{b_{{t - 1 - s},T}^{({t,h})}X_{s,T}}}},} & \left( {2\text{-}33} \right)\end{matrix}$

The next task is to optimize b_(t−1−s,T) ^((t,h)) so that the meansquare prediction error (MSPE), E({circumflex over(X)}_(t+h−1,T)−X_(t+h−1,T))² is minimized. To utilize the properties ofLSW, approximate MSPE can be represented in a wavelet matrix form.)

E({circumflex over (X)} _(t+h−1,T) −X _(t+h−1,T))² ={tilde over (b)}′_(t+h−1,T) B _(t+h−1,T) {tilde over (b)} _(t+h−1),  (2-34)

In Equation 2-34, {tilde over (b)}_(t)({tilde over (b)}_(t−1,T) ⁽¹⁾, . .. , {tilde over (b)}_(0,T) ⁽¹⁾) and B_(1,T) is the covariance matrixconsisting of localized autocovariance whose size is (t+1)×(t+1) and (m,n)-th element is c(n+m/2T, n−m) which is the localized autocovarianceinvolving a wavelet form as mentioned in the background LSW section. Thewhole procedure of forecasting a non-stationary process with long enoughobservations involves two training sets. First, divide the long-termobservations into two training sets and use the first one to generatethe B_(t,T) matrix. Second, use the other training set to optimize{tilde over (b)}_(t) vector such that MSPE is minimized. In practice,the first part involves selection of two more parameters. Once thoseparameters are selected, the model would be ready for prediction.

In Chapter 2, the properties of non-stationary stochastic time seriesand some common methods for modeling or analyzing them are introduced.Forecasting a non-stationary stochastic process can be achieved if theparametric model constructed is proper and robust or the non-parametricmethods are precisely updated.

Chapter 3 Signal Regularity-Based Automated Seizure Prediction Algorithmon Scalp EEG Background and Significance

An epileptic seizure is a transient occurrence of signs and/or symptomsdue to abnormal excessive and synchronous neuronal activity in thebrain. The worldwide prevalence of epilepsy ranges from 0.4% to 1%. Someepileptic patients experience a prodrome or an aura, which can serve asa warning before the signs of seizure onset. Rare patients learn toabort a seizure without external intervention. However, most patientscannot predict or arrest their seizures. In industrialized countries,where antiepileptic drugs and seizure control devices are readilyavailable, about 70% of epilepsy patients are able to gain satisfactorycontrol of their seizures. For patients whose seizures do not respond toantiepileptic medications, less than 50% are candidates for epilepsysurgery. Therefore, approximately 15-20% of epilepsy patients have nochoice but to live their lives with unforeseen and uncontrolled seizureattacks, which cause considerable stress for these patients and theircare-givers and limit the range of daily activities available due tosafety concerns. These lifestyle limitations decrease quality of lifeand may contribute to the increased prevalence of depression in patientswith uncontrolled seizures. If a device could be developed that couldwarn an epilepsy patient of an impending seizure, it could lessen thepsychological stress of epilepsy and improve patient safety.

Numerous studies on intracranial EEG and fMRI data in epilepsy patientssuggest the existence of a preictal transition between an interictal andictal state, and many attempts have been made to identify the mostconsistent patterns of preictal transitions. In the last two decades,signal processing methods in nonlinear dynamics have been popularlyapplied to EEG analyses because of the non-linear nature of neuronalactivity and the paroxysmal character of epileptic seizures. Due to thecomplexity of the implement of prediction, it is not enough to justutilize features showing differences between preictal and ictal states.Decisions about issuing a warning based on features need to be madecontinuously as the EEG signals shift from state to state in the reallife environment. Therefore, many other sophisticated engineeringtechniques derived from statistical learning and control theory werealso employed on EEG data to facilitate seizure prediction performance.Studies utilizing a variety of engineering techniques and also thefeatures on EEG data will be introduced. It is not surprising that mostseizure prediction studies were on intracranial EEG recordings becausethe signals are less influenced by movement, muscle, and other normalphysiologic activities. Chisci, Mavino and coworkers applied classicparametric models and machine learning techniques to accomplish seizureprediction on intracranial EEG from nine patients in the Freiburgdatabase. The authors extracted the parameters of autoregressive modelsfor each EEG channel as dynamic features. The feature values were thenfed into a classifier combining the support vector machine along withKalman filter techniques. Training and testing procedures wereimplemented for 10 runs for each patient. Authors achieved 100%sensitivity and zero false-positive rates on data from two patientswhile the remaining seven patients had false-positive rates among 0.12to 1.03 false warnings per hour. In this study, the definition of afalse positive is when a classifier calls an epoch as preictal and thenclassifies any of the following epochs as interictal before the comingof the next seizure. In contrast, the definition of a true positive ismore intuitive as any epoch classified as preictal before a seizure in asegment containing a seizure. The average time difference betweenwarnings and seizure onsets ranged from 6 to 92 minutes among patients.Netoff, Park and Parhi reported a study that utilized a cost-sensitivesupport vector machine fed with the frequency power of EEG signals toachieve seizure prediction. They applied their method on 219-hour-longEEG data for nine out of the 21 patients in the Freiburg dataset. Theperformance was evaluated using a prediction horizon of five minutes andachieved sensitivity of 77.8% and zero false-positive rates. Mirowski,Madhavan and coworkers applied several bivariate features andclassifiers on the Freiburg database to recognize preictal patterns.Cross-correlation, nonlinear interdependence, dynamical convergence, andwavelet synchrony were estimated in five-second windows. Features werethen aggregated into spatio-temporal patterns to feed into classifiers.Three classifiers were used and compared. The convolutional neuralnetwork outperformed logistic regression and the support vector machine.The best result achieved by using wavelet coherence and a convolutionalneural network was 71% sensitivity with zero false-positive rates for 15out of 21 patients. For those studies mentioned above, they all appliedclassifiers on the Freiburg data and report promising results. It wouldbe intriguing to see if the classifier could work as well on scalp EEGdata. In research conducted by Feldwisch-Drentrup, Schelter, andcolleagues, two dynamic features (mean phase coherence and dynamicsimilarity index) were applied on 1,456 hours of long-term continuousintracranial EEG data from eight patients. The prediction sensitivityincreased (from 25% to 43.2%) when they used logical “AND” combinationsto join the benefits of both features. However, the results of usingindividual methods and logical combination methods were all separatelyoptimized within the same dataset. Therefore, the increase ofsensitivity of combined methods is not surprising. More other Studiesregarding earlier seizure prediction methodologies and development usingintracranial EEG have been extensively reviewed.

Scalp EEG has been used to investigate normal and pathological brainfunction for over 80 years, when Hans Berger published his work andcoined the term “Elektenkephalogram” in 1929. Early attempts to predictseizures from EEG signals began in the 1970s and flourished during the1990s. Because of the non-linear nature of neuronal function and theparoxysmal character of a seizure, non-linear signal processing methodswere popularly applied to EEG for predicting seizures. By applying anon-linear similarity index, Le Van Quyen, Martinerie, and colleagues,studied preictal EEG dynamics in 26 scalp EEG recordings obtained from23 patients with temporal lobe epilepsy. In five patients withsimultaneous scalp and intracranial recordings, changes in thesimilarity index values was observed during preictal period in bothtypes of EEG recordings, and 25 out of the 26 EEG segments showedchanges prior to the occurrence seizures (mean 7 minutes). Although thisstudy did not evaluate the specificity of the similarity index change inlong-term EEG recordings (contained only 50 minutes before seizureonsets), it rendered an encouraging result using only scalp EEG toachieve seizure prediction. Another study by Hively and Protopopescuused L₁-distance and χ² statistic to estimate the dissimilarity indensity functions between the base-windows and the test-windows in 20scalp EEG recordings. The result showed pre-seizure changes in alldatasets with the forewarning times ranged from 10 to 13660 seconds.However, this study only examined one selected channel in each data set,and similar to, specificity of the method was not reported. In twofollow-up studies by the same group, one using all available recordingchannels and the other using a fixed channel, showed similar results.While the results seemed promising, no validation study has beenreported for this method. In the study reported by Corsini, Shoker andcoworkers, 20 sets of simultaneous scalp and intracranial EEG recordingswere analyzed using blind source separation and a short-term Lyapunovexponent (STLmax). This study reported changes of STLmax values beforeseizure onsets and noted that the scalp EEG may give better predictivepower over intracranial EEG when the intracranial electrodes did notrecord the electrical activity in the epileptic focus. However, thepracticality of this method on long-term scalp EEG recordings may belimited due to the lack of an automatic procedure to select the mostrelevant source component. Schad et al. investigated seizure detectionand prediction in 423 hours of long-term simultaneous scalp andintracranial EEG recordings from six epileptic patients. The method usedtechniques based on simulated leaky integrate-and-fire neurons. Thestudy reported that 59%/50% of the 22 seizures were predicted usingscalp/invasive EEGs given a maximum number of 0.15 false predictions perhour. In the study by Bruzzo et al., a small sample of scalp EEGrecordings (115 hours from 3 epileptic patients) was analyzed usingpermutation entropy (PE). By examining the area under the receiveroperating characteristic (ROC) curve, they reported that the decrease ofPE values was correlated with the occurrence of seizures. However, theauthors also concluded that the dependency of PE changes on thevigilance state may restrict its possible application for seizureprediction. More recently, Zandi, Dumont, and colleagues reported aprediction method based on the positive zero-crossing interval series.The method was applied on a 21.5 hour scalp EEG dataset recorded from 4patients with temporal lobe epilepsy. They reported a training result of87.5% sensitivity (16 seizures) with a false prediction rate of 0.28 perhour, where the average prediction time was approximately 25 minutes.James and Gupta analyzed long-term continuous scalp EEG recordings fromnine patients (5 in training set and the other 4 in test dataset). Thedata were processed by a sequence of techniques consisting ofindependent component analysis, phase locking value, neuroscale, andGaussian mixture model. The prediction performance of this methodachieved a sensitivity of 65-100% and specificity of 65-80% as theprediction horizon ranged from 35-65 minutes in the test dataset.

Among these studies, one of the common features used in many of theseizure warning algorithms is the change of synchronization of EEGsignals recorded from different electrode sites. The concept of“synchronization” can be quantified in several different ways such ascoherence, phase synchronization, or entrainment of dynamic features(denoted as “dynamic entrainment” hereafter). In this study, it has beenattempted to identify preictal transitions by detecting dynamicentrainment based on the convergence of PMRS among multiple electrodesites. PMRS is a probabilistic statistic that quantifies the regularityof a time series. It is especially useful when the moment statistics(e.g., mean, variance, etc.) or frequency cannot detect changes in asignal. By further applying a paired t-statistic (denoted as “T-index”hereafter) that quantifies the convergence of two PMRS time series overtime, an automated seizure prediction algorithm may be constructed thatmonitors the change of T-index and issues a warning of an impendingseizure when the T-index curve exhibits the pattern defined by thealgorithm. The general hypothesis is that seizures are preceded by PMRSentrainment; this hypothesis was based upon findings reported previouslyusing a different measure of signal order, STLmax. The predictionparameters and the specific EEG channels to be monitored were determinedby the use of a training dataset. Algorithm performance was thenassessed using an independent test dataset. The performance was furthervalidated by comparison with that from a random warning scheme that didnot use any information from the EEG signals.

Methods Seizure Warning Mechanism Overview

EEG signals were first filtered by a fifth-order Butterworth filter witha band-passing frequency between 1 to 20 Hz (the bandwidth within whichmost ictal epileptiform patterns occur). After the filtering process,for each channel, PMRS was calculated for each non-overlapping5.12-second epoch. Based on the PMRS values, T-indices were thencalculated for each of the selected channel groups. To increase thesensitivity of seizure warning, the proposed algorithm independentlymonitored four T-index curves (i.e., from four channel groups). Awarning was issued when any of the monitored T-index curves metconvergence criteria.

Pattern-Match Regularity Statistic (PMRS)

PMRS is a probabilistic statistic quantifying signal regularity. One ofthe characteristic features of EEG signals during a seizure is therhythmic and regular discharges over a wide range of the brain.Therefore, the first step of the warning algorithm presented in thisstudy calculated PMRS sequentially for each EEG signal analyzed. Therationale of applying this pattern match method (instead of value match)is due to its robustness over scalp signal values, which are usuallymore unstable than their up-and-down trends. The procedure forcalculating PMRS is described below:

Given a time series U={u₁, u₂, . . . , u_(n)} with standard deviation{circumflex over (σ)}{circumflex over (σ_(n))}, a tolerance coefficiente, and a fixed integer m, the two segments in U (x_(i)={u_(i), u_(i+1),. . . , x_(i+m−1)}, x_(j)={u_(j), u_(j+1), . . . , u_(j+m−1)}) areconsidered pattern-matched to each other when Equation 3-1 is fulfilled.

{|u _(i) −u _(j) |≦e{circumflex over (σ)}{circumflex over (σ_(n))}}

{|u_(i+m−1) −u _(j+m−1) |≦e{circumflex over (σ)}{circumflex over(σ_(n))}}

{sign(u _(i+k) −u _(i+k−1))=sign(u _(j+k) −u _(j+k−1)),k=1,2, . . .,m−1}  (3-1)

In Equation 3-1, the first two criteria require value match to someextent at both the beginning and ending points of two segments, where ewas set to be 0.2 empirically. The third criterion requires patternmatch between x_(i) and x_(j) within a range of m (set as 3 in thisstudy). To calculate PMRS, first define a conditional probability,p_(i).

p _(i) =Pr{sign(u _(i+m) −u _(i+m−1))=sign(u _(j+m) −u _(j+m−1))|x _(i)and x _(j) are pattern match},  (3-2)

Given m, p_(i) can be estimated as {circumflex over (p)}_(l) as inEquation 3-3.

$\begin{matrix}{{\hat{p_{l}} = \frac{\# \mspace{14mu} {of}\mspace{20mu} \begin{Bmatrix}{\left\lbrack {x_{j}^{\prime}s\mspace{14mu} {pattern}\mspace{14mu} {match}\mspace{14mu} {with}\mspace{14mu} x_{i}} \right\rbrack\bigwedge} \\\left. \left\lbrack {{{sign}\left( {u_{i + m} - u_{i + m - 1}} \right)} = {{sign}\left( {u_{j + m} - u_{j + m - 1}} \right)}} \right) \right\rbrack\end{Bmatrix}}{\# \mspace{14mu} {of}\mspace{14mu} \left\{ {x_{j}^{\prime}s\mspace{14mu} {pattern}\mspace{14mu} {match}\mspace{14mu} {with}\mspace{14mu} x_{i}} \right\}}},\mspace{20mu} {1 \leq j \leq {n - m}}} & \left( {3\text{-}3} \right)\end{matrix}$

In Equation 3-3, 1≦i≦n−m. Finally a PMRS can be estimated.

$\begin{matrix}{{PMRS} = {{- \frac{1}{n - m}}{\sum\limits_{i = 1}^{n - m}{\ln \left( \hat{p_{l}} \right)}}}} & \left( {3\text{-}4} \right)\end{matrix}$

As the time series U develops into a more regular state, {circumflexover (p)}_(l)s become larger and PMRS decreases as a result.

PMRS Convergence (T-Index)

T-index is basically the paired t-statistic function used to quantifythe degree of convergence between two PMRS time series. For two timeseries X_(i) and X_(j) (the PMRS value time series), their values in acalculation window W^(t) with a size of n data points are presented asL_(i) ^(t) and L_(j) ^(t) in Equation 3-5 and Equation 3-6.

L _(i) ^(t) ={X _(i) ^(t) ,X _(i) ^(t+1) , . . . ,X _(i)^(t+n−1)}  (3-5)

L _(j) ^(t) ={X _(j) ^(t) ,X _(j) ^(t+1) , . . . ,X _(j)^(t+n−1)}  (3-6)

Then the pair-wise differences between L_(i) ^(t) and L_(j) ^(t) can beexpressed in Equation 3-7.

D _(ij) ^(t) =L _(i) ^(t) −L _(j) ^(t) ={X _(i) ^(t) −X _(j) ^(t) ,X_(i) ^(t+1) −X _(j) ^(t+1) , . . . X _(i) ^(t+n−1) X _(j) ^(t+n−1) }={d_(ij) ^(t) ,d _(ij) ^(t+1) , . . . ,d _(ij) ^(t+n−1)},  (3-7)

The T-index over the calculation window W^(t) between the two timeseries is calculated using Equation 3-8.

$\begin{matrix}{{Tind}_{ij}^{t} = \frac{{\overset{\_}{D}}_{ij}^{t}}{{\hat{\sigma}}_{D_{ij}^{t}}/\sqrt{n}}} & \left( {3\text{-}8} \right)\end{matrix}$

In Equation 3-8, D _(ij) ^(t) and {circumflex over (σ)}_(D) _(ij) _(t)are the sample mean and the sample standard deviation of D_(ij) ^(t),respectively.

PMRS convergence is a process during which one EEG signal is influencedby or coupled with another with respect to the signal regularity. Thisphenomenon was used as a dynamical pattern for warning of an impendingseizure. FIG. 3-1 shows the PMRS traces derived from three EEG channels(F8, T4, and T6) (upper panel) and their average T-index values (bottompanel) over a 350-minute interval containing a seizure. As shown in thePMRS plot, all PMRS values of the three channels dropped seconds afterthe seizure onset (indicated by a black vertical line), which was due tothe extreme signal singularity during the ictal period. Moreimportantly, the PMRS values became convergent about 60 minutes beforethe seizure onset, which caused the decrease of T-index values (shown inFIG. 3-1).

Selection of Channel Groups for Seizure Warning Monitoring

Channel groups were selected such that they were bilaterally symmetricalong the midsagittal line. In order to avoid frequent eye movementartifact, electrodes Fp1 and Fp2 were excluded. In addition, to avoidthe regular alpha rhythm pattern that might give potential falsepositive warnings because of the nature of the regularity statisticsused in the prediction algorithm, electrodes O1 and O2 were alsoexcluded. The frontal electrodes F3 and F4 were included in theanalysis, whereas electrodes F7 and F8 were not included in order tominimize muscle artifact due to chewing.

The four channel groups selected in this study resemble the well-knownand popular anterior-posterior bipolar (“double banana”) montage. Thesefour channel groups were: (F7, T3, T5), (F3, C3, P3), (F4, C4, P4), and(F8, T4, T6). The main rationale for this selection was that mostseizures occurring in patients with temporal lobe epilepsy start from achannel group on one side (hemisphere) of the brain, and therefore,intuitively, these channels were more likely to be entrained with eachother during the preictal transition period. By monitoring four channelgroups that cover both hemispheres, the algorithm was enabled to detecta PMRS convergence that preceded an impending seizure initiated fromeither hemisphere.

For each channel group k, k=1-4, the group T-index over the calculationwindow W^(t) is denoted as:

GTind_(k) ^(t)={Σ_(i=1) ²Σ_(i′=i+1) ³ Tind_(ii′) ^(t)}/3  (3-9)

The window size for calculating one group T-index is 60 PMRS datapoints, with 59 points overlapping from t to t+1 (i.e., sliding window).The warning algorithm only monitored the four group T-index curvesinstead of individual pair-wise T-indices.

Detection of PMRS Convergence

The proposed prediction algorithm, instead of attempting to detect asignal threshold crossing, was set to detect a certain pattern ofT-index dynamics that gradually descends from a baseline value.Therefore, the first step to detect such a pattern was to determine anupper threshold U_(T) (i.e., baseline value, see an illustration in FIG.3-2). A decrease of a group T-index from U_(T) was considered as adesired condition for a potential PMRS convergence, and a convergencewas identified when the group T-index values fell below a lowerthreshold L_(T). It was further defined that a period of convergenceshould be maintained for at least several minutes. For each groupT-index, its U_(T) was set to be the asymptotic 95^(th) percentile inthe preceding 12 minutes, as described in equation 10. The duration of12 minutes was decided by training across multiple patients in thetraining dataset.

U _(T)=mean(T _(12min)(t))+2×std(T _(12min)(t)),  (3-10)

However, if half of the following 16 T-index values (representing 81.92seconds of data) were greater than U_(T), U_(T) was updated as themedian of these 16 values. This operation is described in equations 11and 12.

$\begin{matrix}{{U_{T} = {{U_{T} \times \left( {1 - {I(A)}} \right)} + {{median}\; (T) \times {I(A)}}}},} & \left( {3\text{-}11} \right) \\{{A = \left\{ {{\sum\limits_{k = 1}^{16}{H\left\lbrack {{T\left( {t + k} \right)} - U_{T}} \right\rbrack}} > \frac{16}{2}} \right\}},} & \left( {3\text{-}12} \right)\end{matrix}$

In Equation 3-12, I(.) is an indicator function, A is the criterion ofupdating U_(T), and T(t) is the group T-index value at time t. If A isfalse, then U_(T) is kept as it is. H(.) denotes the Heaviside stepfunction. U_(T) is not updated with a maximum value because theexistence of artifact in raw EEG recordings could affect the T-indexcausing an abrupt surge and therefore produce an abnormally high U_(T).Once U_(T) is determined, L_(T) is equal to U_(T) minus D, where D is aparameter that would determine the sensitivity of the algorithm. Ingeneral, the bigger D is, the less sensitive the algorithm will be, butthe less susceptible the algorithm will be to false warnings.

In addition to the above detection rules, the period of a descendingT-index pattern from U_(T) to L_(T) should continue for more than t_(t)minutes to be regarded as a PMRS convergence. Therefore, the algorithmwould issue a warning of an impending seizure only when any of themonitored T-index curves traveled from U_(T) to L_(T) with the travelingtime longer than t_(t) to ensure a gradual descendent pattern.

Once a convergence was identified, other convergences that followed(identified within the seizure warning horizon (SWH) by the warningalgorithm using the same parameter settings) were silenced due to thepossibility that a convergence preceding an onset may be much longerthan the t_(t) value and cause several convergence identifications. Thework flow is depicted in FIG. 3-3 and gives an overview of the seizurewarning algorithm.

EEG Data Characteristics Subjects and EEG Recording Specifications

All subjects were 18 years of age or older admitted to either AlleghenyGeneral Hospital (AGH, Pittsburgh, Pa.) or the Medical University ofSouth Carolina (MUSC, Charleston, S.C.) for inpatient seizure monitoringfor diagnostic purposes or presurgical evaluation. Data collectionprocedure was approved by the Investigational Review Boards of AGH andMUSC, and the Western Investigational Review Board (WIRB). The EEGrecordings at MUSC were obtained using XLTEK monitoring systems(Oakville, Ontario, Canada) with a sampling rate of 256 Hz and the EEGrecordings at AGH used 128-channel Nicolet BMSI-6000 systems (Viasys,Madison, Wis., USA) with a 400 Hz sampling rate. The EEGs recorded atboth institutions used a referential montage and the 19-electrodeinternational 10-20 system of electrode placement. The exact locationsof referential electrodes placed in the dataset were decided on-site andusually followed the recommended location of the Cz and Pz electrodes assuggested by the American Clinical Neurophysiology Society. All segmentswere reviewed by the collaborative clinical sites (AGH and MUSC) and allseizure events were verified by KK and JH, respectively.

Data Selection

Because a scalp EEG might be severely contaminated with artifact,typical muscle contraction, and movement artifacts include blinking,chewing, and talking, only EEG segments with a tolerable level ofartifact, i.e., not present for more than 50% of the recording and notinvolving more than 50% of the recording channels, were included in thisstudy. There were a total of 71 EEG recording segments from 71 patientsselected for this study. All EEG segments were long-term recordings(mean=25.16 hours) containing at least one seizure. EEG segmentscontaining more than one seizure within any two-hour interval wereexcluded to avoid potential overlap between preictal and postictalperiods of consecutive seizures. Thus, the warning algorithm had asufficiently long period of observation to detect the transition frominterictal to ictal states. The dataset was randomly divided intocomplementary training (n=35) and test (n=36) sets. The total EEGrecordings used in this study was 1786 hours and contained 98 seizures.The individual recording durations of each subject in the whole dataset(both training dataset and test dataset) are shown in FIG. 3-4.

Statistical Evaluation Estimation of Performance Statistics

The performance of the seizure warning algorithm varies with the lengthof the seizure warning horizon (SWH)—the longer the SWH, the lower thesensitivity and false positive rate (FPR). In this work, SWH was definedas the time window following a seizure warning during which the patientwas likely to have a seizure. It is worth noting that Winterhalder etal. defined the SWH (same as seizure prediction horizon, SPH)differently as a short intervention preparation period. The definitionof SWH here is closer to the “seizure occurrence period” defined by thisgroup. A seizure warning was considered true only when at least oneseizure onset occurred within the SWH following the warning. Otherwise,the warning was considered a false positive.

After determining the outcome (i.e., true or false) of each warning,sensitivity was estimated by dividing the number of correct warnings bythe total number of seizures in the EEG segments, i.e. the proportion ofseizures correctly predicted (within SWH). FPR (per hour) was calculatedby dividing the total number of false positives by the total recordinghours outside the SWH before each seizure. The SWH period before eachseizure was excluded because it was impossible for a false warning tooccur within the SHW before each seizure.

Performance Statistical Validation—Comparison with a Random SeizureWarning Scheme

Several methods have been proposed for validating the performance of aseizure warning algorithm. Each study tested a specific statisticalhypothesis. In this study, the prediction sensitivity estimates werecompared with a random prediction scheme that issued random seizurewarnings (in time) during the recording. Under this random scheme, theinterval between two consecutive warnings followed an exponentialdistribution, with a condition that the length of any interval could notbe smaller than the SWH. This additional condition ensured that the twocompared algorithms (test algorithm and prediction scheme) wereconsistent in terms of the management of nearby warnings (within theSWH). Although this condition diminished some randomness from the randomprediction scheme, it provided a meaningful comparison of predictionperformance between the two prediction methods.

In order to have an unbiased comparison for sensitivity, two parameterswere controlled: 1) length of the SWH, and 2) total number of warnings.For each of the test EEG segments, the total number of warnings allowedby the random predictor was set to be the same as that issued by thetest algorithm. For example, if the test algorithm issued two truepositive and two false warnings in a specific segment, the randompredictor would only randomly issue a total of four warnings. Under thesame number of warnings allowed, the study compared overall predictionsensitivity across all test patients. A distribution of overallsensitivity by the random prediction scheme was generated from 1000simulations, and was used to assess how significant the predictionperformance by the test algorithm was over the random scheme.

It is worth mentioning that imposing the above conditions in the randomwarning algorithm compared in this study allowed for comparison with theproposed algorithm. Since this was a conditioned random process, thespecific null hypothesis tested in this study was: “The overallsensitivity (or FPR) of the test algorithm is the same as that of arandom prediction scheme with the same number of warnings issued by thetest algorithm.”

Cross-Validation on Performance Comparison

The comparison of performance between the proposed warning algorithm andthe random scheme was further implemented through cross-validation sothat the comparison result is more independent from the constituents inthe training dataset and test dataset. To execute the cross-validate,all 71 segments may be randomly divided into training and test datasetsfor 20 times to generate 20 individual trials and then repeated thetraining and test procedures on each trial. In each trial, theparameters (t_(t) and D) were selected as those causing the performanceto exceed 70% sensitivity and resulting in the least FPR in the trainingdataset. The parameters were input to the test algorithm run on testdataset to generate the final performance results in the trial. For eachtrial, the final performance result of the test dataset was comparedwith the random seizure warning scheme mentioned above. Please note thateach test segment in each trial can have a different number of warningsand even the same segment in different trial would have a differentnumber of warnings if the parameters used in the two trials weredifferent. Using fixing parameters for a group of epileptic patientswith different seizure types and physiological backgrounds can hardlyoutperform that using best setting parameters for each individualpatient. However, these results can offer a perspective on therobustness of using the test algorithm among several patients when theprevious individual patient data is not available or sufficient fortraining.

Combining p-Values

In each trial, the comparison result between the test algorithm and therandom scheme is quantified as a p-value, and indicates how significantthe prediction performance by the test algorithm was over the randomscheme. After obtaining all the p-values from the 20 cross-validationtrials, the overall statistical significance of comparison between thetwo seizure prediction schemes was estimated by mapping the p-values toa z distribution. The z-transform method was chosen for testing combinedp-values because it is not differentially sensitive to data that supportor refute a common null hypothesis, whereas Fisher's test, anotherwidely-used method for testing combined p-values, is more sensitive todata refuting a common null hypothesis. All p-values of each trial weretransformed to a z-value, z_(p), in a N(0.1) distribution such thatPr[N(0,1)≦z_(p)]=p-value. The overall z-value, Z_(overall), is carriedout by calculating Σz_(p)/√{square root over (k)}. Lastly, Z_(overall)is transformed back to a combined p-value.

Results Training Results

The proposed algorithm was applied on the training dataset (n=35) tooptimize settings for the parameters D and t_(t). Parameter D requeststhe minimal decrease that a group T-index value needs to fall from U_(T)to be considered as a PMRS convergence, and t_(t) gives the constrainton the minimal time length for a group T-index curve traveling fromU_(T) to L_(T). One of the training results is presented in FIG. 3-5,which shows the receiver operating characteristic (ROC) curves(sensitivity vs. false positive rate) under different t_(t), rangingfrom 10 to 30 minutes with 10-minute increments. Each ROC curve containseight settings of parameter D, each representing the performance under acertain D value, ranging from 1 to 8 (in T-index units) with incrementsof 1. As D increased, each ROC curve first reached its peak sensitivityand then declined as D kept increasing. This was due to the constraintof t_(t): when D was set to a small value, T_(t), the time intervalduring which a group T-index curve drops from U_(T) to L_(T), was likelyto be shorter than the minimally required t_(t), and therefore the dropof the group T-index curve could not be considered as a PMRSconvergence. As D became bigger, the sensitivity started increasinguntil the D value became too large to start affecting the number ofgroup T-index drops (i.e., sensitivity started to reduce).

For the training processes of all trials, the best performance wasdefined as when the proposed algorithm achieves sensitivity of at least70% with the lowest possible false positive rate. The best parameterconfiguration was optimized within each trial according to thedefinition. As a result, the optimized prediction parameter values weredecided for each trail after reviewing all performances using a fullrange of possible parameter settings. These optimized parameters werefixed for the performance evaluation in the paired test dataset withinthe same trail.

Test Results

The test dataset consisted of 36 long-term EEG segments from 36 patientsin each trial. With the parameter settings determined in the trainingdataset of the same trial, the test algorithm gave sensitivities among0.6 and 0.72 (mean=0.65) while the false-positive rate was among 0.22and 0.26 (mean=0.24).

The random warning scheme was designed to issue random warnings with thesame number as that issued by the test algorithm for each segment in thetest dataset. After the random scheme ran through all of the patients inthe test dataset an overall sensitivity was calculated. The randomscheme repeated 1000 times on the test dataset and thus generated 1000overall sensitivity estimates. The overall sensitivities of all 20trials are listed in the Table 3-1. The overall sensitivity achieved bythe test algorithm was then compared with the sensitivity distributionof the random warning scheme. For example, the performance comparison inone of the trials is shown in FIG. 3-6. The combined p-value over 20cross-validation trials showed that the proposed warning algorithmachieved a better performance (p=0.015) than the compared random warningscheme.

Discussion Dataset Requirement for Performance Evaluation

Amongst the early seizure prediction studies, short EEG recordingcontaining seizure onsets were often used. The claim of sensitivity canbe very satisfactory especially when the result was optimized for allthe data without separating them into training and test dataset.Although these studies were absolutely valuable as pioneers in thefield, simply interpreting those results as successfully predictingseizures would be too optimistic. Concerns of overestimation of resultssoon developed and drew much attention amongst many research groups.Using interictal EEG recordings and a separate test dataset forperformance evaluation was advocated to fairly estimate the falsepositive rate and avoid misleading conclusions. In sum, potentialmisleading factors that could overrate the performance should be ruledout. These concerns should not be regarded as unduly skeptical orcritical of the rigor of these early study designs; rather, theyrepresent the positive intention of expanding the previous promisingresults to more challenging test scenarios closer to reality. In thisstudy, for the interictal recording requirement, 1786-hour-long EEGsegments in total were included, which included 98 seizures in 71segments. Each seizure is at least three hours away from each other. Ifone defines a preictal period as three hours before the onset of aseizure, the dataset contains 1492 hours of interictal period. Theseinterictal EEG signals should be sufficient for a reasonable estimate ofthe false-positive rate. For the separate testing dataset requirement,the database was randomly divided into training and test dataset for 20trials. Every EEG segment was recorded from a different patient andbelongs either to the training or test dataset in a trial. This schemeof dataset treatment avoids both in-sample and patient-specificoptimization.

Patient-Specific Optimization

If one has a long enough EEG recording, containing several seizures foreach patient in the dataset, patient-specific training and a test schemecan be implemented, and sensitivity as well as specificity can still befairly estimated. For many studies using classifiers to predictseizures, the classifier is usually optimized using a portion of the EEGsegment of one patient and then tested on the rest of EEG segments fromthe same patient. Training parameters vary from method to method. Aslong as all parameters were fixed once after the training procedureends, no more optimization or changing of parameter values should beallowed in the test process; otherwise, the result should not be viewedas an independent test result but a result of further in-sampleoptimized procedure. Although the dataset is divided into two parts, thefurther optimized result should not be viewed as an independent testresult since some parameters have fixed before the further optimization.In a recent study by Feldwisch-Drentrup et al., they divided theirdataset into two sets. They executed channel selection in the first partof the data, and optimized prediction intervention time and thresholdsfor two individual methods, as well as two logical combinations of twomethods in the remaining part of the dataset. The results of the fourprediction methods were all in-sample optimized. In the study, thealgorithm was tested using the same parameters across many patients forone trial and predefined a fixed prediction horizon for all patients.These test results are valid as long as the test datasets have longenough intracranial segments for estimating a false-positive rate andenough seizure onsets for estimating sensitivity. The key point is tointerpret the results with the design and the hypothesis of the studycombined. The study is intended to test the algorithm across patientsagainst a random predictor issuing the same number of warningsuniformly.

Vigilance States as a Possible Confounding Factor

In addition to using interictal segments and separate test datasets forperformance evaluation, there are other concerns related to studydesign. One possible issue is the control of other confounding factors.For example, the scalp EEG signals respond to vigilance states moreobviously than intracranial EEG signals. Chewing and talking muscleartifacts on the temporal region electrodes during awake state, andprevailing alpha waves on occipital region electrodes during eyes-closedbut awake state, are very dominant in the scalp EEG signals. If theinterictal segments were all recorded during the awake state whileseizures occurred during sleep, an algorithm detecting drowsy state canachieve a high sensitivity predicting seizures. One way to rule out ofthis erroneous relationship is selecting many interictal and preictalrecordings under the same vigilance state for both training and test.Whether the algorithm is more sensitive to a vigilance state or anepileptic state change can be further clarified by conducting anothercontinuous study. However, the existence of many different vigilancestates in the data allows one to test the robustness of a proposedalgorithm under many vigilance states. In this study, all of the EEGsegments are continuous long-term EEG segments and most of them are longenough to go through several vigilance states. A seizure could happenduring the sleep or awake states in a continuous recording. By includingmany continuous long-term segments from many different patients,seizures were scattered among several vigilance states instead of oneprevailing condition.

Significance Test

The result of this study should not be interpreted to mean that theperformance of the proposed algorithm outperformed any random predictionscheme. This study tested a specific null hypothesis based on the randomprediction scheme designed. In fact, a comparison of performance with acompletely random prediction scheme may have very little meaning, bothscientifically and practically. Certain conditions should be imposed onthe random scheme to prevent potential bias in the comparison,especially for the assessment of a sophisticated EEG-based predictionalgorithm. Nevertheless, there are several studies that have applieddifferent random warning schemes or surrogate data to test differenthypotheses. The comparison of the proposed algorithm with other randomwarning schemes should be followed. Furthermore, studies on a largersample size and the assessment of the proposed algorithm with differentSWHs should be conducted in the future.

TABLE 3-1 The sensitivities, false-positive rates and p-values achievedby the proposed algorithm on training and test dataset. Thesensitivities and FPRs of the training datasets presented in this tableare the best performances, defined as reaching sensitivity over 70%,meanwhile, resulted in the least FPR. The sensitivities and FPRs of thetest datasets resulted from applying the proposed algorithm using thebest parameter configuration found using the training dataset. Thep-values in the test datasets were estimated by comparing the testalgorithm to the random warning scheme (1000 simulations in each trial)that issued the same number of warnings. Training Dataset Test DatasetTrial Number Sensitivity FPR (times/hour) Sensitivity FPR (times/hour) 10.71 0.22 0.62 0.24 2 0.71 0.27 0.66 0.26 3 0.73 0.25 0.7 0.24 4 0.720.23 0.71 0.25 5 0.76 0.25 0.62 0.23 6 0.7 0.23 0.62 0.24 7 0.74 0.280.67 0.26 8 0.71 0.24 0.61 0.22 9 0.72 0.24 0.6 0.23 10 0.76 0.25 0.650.25 11 0.76 0.24 0.67 0.25 12 0.75 0.24 0.57 0.22 13 0.77 0.24 0.680.22 14 0.71 0.24 0.6 0.23 15 0.7 0.24 0.72 0.25 16 0.72 0.24 0.65 0.2617 0.7 0.26 0.7 0.24 18 0.72 0.23 0.61 0.24 19 0.72 0.24 0.67 0.23 200.72 0.24 0.7 0.25 Average 0.73 0.24 0.65 0.24

Chapter 4 Psycogenic Non-Epileptic Seizure and Complex Partial SeizurePatients Classification Preface

About 25% of the patients visiting an epilepsy monitoring unit arediagnosed as psychogenic nonepileptic seizure (PNES) patients and needto be distinguished from other patients having epileptic seizures inorder to be correctly treated. Patients have PNES onsets that can beidentified by monitoring the ictal pattern shown either in scalpelectroencephalogram (EEG) or video. The classification using ictalsymptoms requires several seizing onsets from a patient. The number orfrequency of seizure onsets during a patient's stay in EMU cannot beprecisely anticipated in advance. It would be convenient for EMUschedule and cost efficiency if one could identify PNES patients usingonly awake and interictal scalp EEG recordings. This study researchedthe connectivity and features of awake and relaxed interictal EEGsignals. The subjects include seven patients having PNES and anotherseven patients having complex partial seizures (CPS) with fixed foci.

The hypothesis is that the inter-hemispheric connectivity and networkdynamic features from the EEG of a CPS patient are different from thatof a PNES patient because the repetitive focal discharges impair morearound the epileptogenic zone rather than the contralateral hemisphere.Conversely, patients with PNES onsets should not have this abnormalphysioelectrical consequence, and are supposed to have a sounderfunctional network or more similar values of features between bothhemispheres.

Small-World Network in EEG Data

A brain is a structurally and functionally complex network of neurons.The functional network reflects the connectedness among brain regions interms of neuronal activity. Graph theoretical analysis is a mathematicaltool to reveal topological characteristics of a network. Applying graphtheoretical analysis on the EEG data reveals the brain functionalnetwork features. One of the network structures, called “small-worldnetwork” can be identified, when there is a balance between localindependence and global integration in the network. The balance can beevaluated by quantifying two graph features, called the local clusteringcoefficient and the characteristic path length. A small world networkhas a relatively high cluster coefficient and a small characteristicpath length. Small world networks are usually compared to a network witha lattice-like configuration or to a random network; the two extremecases. A regular lattice network is characterized by a relatively highcluster coefficient and a long average path length. On the other hand, arandom network has a relatively lower cluster coefficient and a shorteraverage path length. Small-world networks are efficient at informationprocessing, cost-effective, and are relatively resilient to networkdamage and, as a result, can be regarded as the ideal model for anormally functioning brain network.

Graph Theoretical Analysis on EEG Signals Preprocess

While reviewing the EEG data, one can easily observe that EEG data hasone salient characteristic, the rhythmic oscillation. Not surprisingly,therefore, frequency analysis has been broadly used to preprocess andanalyze EEG signals. It is widely accepted that the different strengthsof the frequency components can reveal different brain states.

In analyzing EEG signals, filtering is usually implemented to removeartifacts or to help focus on the frequency bands of interest. Beforedoing a network analysis of EEG data it is important to remove theartifacts that contribute to the multiplication of channels. Withoutdoing this first, the shared input from the artifacts would likelyresult in the existence of a spurious association amongst channelshaving the same artifact. Therefore, applying a power line frequencynotch filter on EEG data is highly recommended. Any other artifact thatwould affect more than one channel, such as movement of the referenceelectrode, should be removed; otherwise, the result of graph should becarefully inspected to rule out the spurious coupling strength amongstvertices with common artifact sources. Some montages can be used toeliminate possible common artifacts between two channels. Beside usingmontage, average-reference is also commonly used when dealing with EEGdata.

In addition to using conventional fast Fourier transform-based filteringtechniques to obtain EEG signals within a range of frequencies, somestudies implemented wavelet transform to extract components with morespecific time-frequency resolution. For example, Palva et al. (2010)used Morlet wavelets to filter EEG data into 36 frequency bands beforeconstruct the oscillatory phase synchronized network from EEG data.

Revealing Network Structures in EEG Data

The brain functional network can be represented as a graph. Graphtheoretical analysis is then applied after the graph has beenestablished. To do so, first define vertices and edges in the EEG data.If the EEG channels are designated as the vertices of a graph, an edgebetween two vertices signifies a functional connection between two EEGchannels. One would expect a larger correlation between two EEG channelswhen there is an edge between them. Edges can also be values quantifyinghow well the two vertices correlate in weighted graphs. Most of thestudies using graph theoretic analysis on EEG data assume thatstatistical interdependencies between EEG time series reflect functionalinteractions between neurons in the brain regions. There are manystatistical metrics computing the degree of association between timeseries. Some of the most commonly used metrics in EEG functionalconnectivity graph studies will be discussed below.

Phase locking value (PLI) is a statistic measuring thefrequency-specific synchronization between two neuroelectric signals.This is a method focusing on the phase information of time series and isdifferent from coherence, which gives the interrelation of bothamplitude and phase between signals. The PLI is calculated by firstextracting the instantaneous phase information of signals through eitherwavelet transform or Hilbert transform. Both methods lead to a similarresult in real world EEG data. When doing the wavelet transform on asignal x(t), a wavelet function, ψ(t), is first chosen. Gabor and Morletwavelet functions have both been applied on EEG data. Then the waveletcoefficient time series, W_(x)(t) can be computed by convolute x(t) andW_(x)(t).

W _(x)(t)=∫Ψ(t′)x(t′−t)dt′  (4-1)

Then the phase time series, φ(t), can be computed.

$\begin{matrix}{{\varphi^{Xi}\left( {n,s} \right)} = {\arctan \frac{{imag}\left( {W^{Xi}\left( {n,s} \right)} \right)}{{real}\left( {W^{Xi}\left( {n,s} \right)} \right)}}} & \left( {4\text{-}2} \right)\end{matrix}$

When doing the Hilbert transform, the phase information, φ(t), of asignal s(t) is obtained through Equation 4-3.

$\begin{matrix}{{\varphi (t)} = {\arctan \frac{\overset{\sim}{s}(t)}{s(t)}}} & \left( {4\text{-}3} \right) \\{{\overset{\sim}{s}(t)} = {\frac{1}{\pi}{p.v.{\int_{- \infty}^{\infty}{\frac{s(\tau)}{t - \tau}\ {\tau}}}}}} & \left( {4\text{-}4} \right)\end{matrix}$

In Equation 4-4, p. v. denotes the Cauchy principal value in theequation. Then the relative phase, φ_(1,1)(t), can be calculated asEquation 4-5.

φ_(e,ƒ)(t)=φ_(e)(t)−φ_(ƒ)(t)  (4-5)

The subscript n and m signify the relative phase relationship betweenchannel e and f. Finally, it is possible to compute the PLI, PLV, usingEquation 4-6.

$\begin{matrix}{{PLV} = {{\frac{1}{N}{\sum\limits_{j = 1}^{N}\; ^{_{\phi_{,f}}{({{j\Delta}\; t})}}}}}} & \left( {4\text{-}6} \right)\end{matrix}$

There are other computational details including widowing and unwrappingthe instantaneous phase being considered. Interested readers can referto those references mentioned above. PLI was used in a study byDimitriadis et al. to construct the functional connectivity graphs from30-electrode EEG data. They utilized surrogate data to generate abaseline distribution of random PLIs and then determined the functionalconnections (edges) if there was a significantly different (p<0.001) PLIfor a specific pair. The surrogate data is generated by permuting theorder of trials of one signal repeatedly.

Synchronization likelihood (SL) is a statistic measuring the non-linearsimilarity between time series. SL offers an extended perspective ofcorrelation that is not limited in terms of linear relationship; whilecoherence, has the limitation of rendering only the linear correlationas a function of frequency. The computation of SL first requires theinterested time series, X_(k,i) (k=1, . . . , M as channel number andi=1, . . . , N as time indices), to be reconstructed using the method ofdelays.

X _(k,i)=(x _(k,i) ,x _(k,i+l) ,x _(k,i+2l) , . . . ,x_(k,i+(m−1)l))  (4-7)

In Equation 4-7, 1 denotes the lag and m is the embedding dimension. Anumber H_(i,j) is then defined to denote the number of channels X_(k,i)and X_(k,j) that are closer than a crucial distance, ε_(k,i).

$\begin{matrix}{H_{i,j} = {\sum\limits_{k = 1}^{M}\; {\theta \left( {ɛ_{k,i} - {{X_{k,i} - X_{k,j}}}} \right)}}} & \left( {4\text{-}8} \right) \\\left\{ \begin{matrix}{{{{{if}\mspace{14mu} {{X_{k,i} - X_{k,j}}}} \geq {ɛ_{k,i}:S_{k,i,j}}} = 0}\;} \\{{{{if}\mspace{14mu} {{X_{k,i} - X_{k,j}}}} < {ɛ_{k,i}:S_{k,i,j}}} = \frac{H_{i,j} - 1}{M - 1}}\end{matrix} \right. & \left( {4\text{-}9} \right) \\{S_{k,i} = {\frac{1}{2\; N}{\sum\limits_{j = 1}^{N}\; S_{k,i,j}}}} & \left( {4\text{-}10} \right)\end{matrix}$

We can see that SL, S_(k,i), is a measure describing how well channel kat time i is correlated to all other M−1 channels. Ponten et al. usedSynchronization Likelihood (SL) as the connectivity metric in theirstudy about the relationship between epilepsy and small-world networks.The edges between vertices were determined by applying a threshold onthe SL value between two time series recorded at the two vertices withinan epoch. In their study, for the threshold, it was decided to keep thegraph density as a constant. In other words, the threshold increasesslightly until the average number of edges of each vertex in a graph isequal to a certain number.

Non-linear Independence Measure (NIM) estimates the strength offunctional coupling of two time series embedded in their respectivephase-spaces. Let Γ_(n,j) and Λ_(n,j), j=1, . . . , w, denote the timeindices of the w nearest neighbors of two reconstructed vectors, X_(n)and in a phase space respectively. For each X_(k,i), the mean squaredEuclidean distance to its w neighbors can be defined as Equation 4-11.

$\begin{matrix}{{R_{i}^{(w)}(X)} = {\frac{1}{w}{\sum\limits_{j = 1}^{w}\; \left( {X_{k,i} - X_{k,\Gamma_{i,j}}} \right)^{2}}}} & \left( {4\text{-}11} \right)\end{matrix}$

We can also define Y-conditioned mean squared Euclidean distance asEquation 4-12 but replace the nearest neighbors by the equal timepartners of the closest neighbors of X_(k′,i).

$\begin{matrix}{{R_{i}^{(w)}\left( {XY} \right)} = {\frac{1}{w}{\sum\limits_{j = 1}^{w}\left( {X_{k,i} - X_{k,\Lambda_{i,j}}} \right)^{2}}}} & \left( {4\text{-}12} \right)\end{matrix}$

The NIM, N(X|Y), is computed as Equation 4-13.

$\begin{matrix}{{N^{(w)}\left( {XY} \right)} = {\frac{1}{N}{\sum\limits_{i = 1}^{N}\; \frac{{R_{i}(X)} - {R_{i}^{(w)}\left( {XY} \right)}}{R_{i}(X)}}}} & \left( {4\text{-}13} \right)\end{matrix}$

The asymmetry of NIM (N(X|Y)˜=N(Y|X)) is the main advantage of thisnonlinear measure because it can deliver directional information betweenvertices. In the study done by Dimitriadis et al., they used both SL andNIM to quantify the extent of association between vertices. In theirstudy, they found that the correlation dimension was around six from thesleeping EEG data of 10 healthy subjects. Furthermore, SL and NIM areweakly correlated and considered as revealing complementary informationregarding the functional connectivity.

Phase lag index (PLAI) quantifies the asymmetry of the distribution ofinstantaneous phase differences between two time series. If the relativephase time series between channel e and f, φ_(e,f)(t) have beencalculated, then the PLAI can be defined as |

sign[φ_(e,ƒ)(t)])|

|(

.

signifies the expectation value operator). PLAI ranges from 0 to 1. PLAIvalues greater than 0 suggest the existence of phase locking to someextent, and values equal to 0 signify no coupling or no coupling with aphase difference centered around 0±π radians. PLAI is supposed to ruleout the synchronization due to instantaneous volume conduction or acommon source that is the main cause giving spurious synchronization.

Results of Studies Using Small-World Network Analysis to EEG Data

The existence of small-world structure in the cortex of a human brainhas been observed through many measuring methods including fMRI,magnetoencephalography (MEG), and EEG. Many studies use graphtheoretical analysis to distinguish pathological consequences oridentify the state of human brains. The following paragraphs introducesummaries of some interesting studies related to small word networks inEEG data. Readers interested in more detail about the application ofsmall world networks in EEG will find more information following thereferences.

An ageing study found out that the clustering coefficient and the valueof path length were both lower in the elderly compared to the youngsubject group. This implies that the brains of the elderly subject groupappear to be closer to random networks.

Small-world networks are supposed to be very efficient for datatransfer. Epilepsy as a disease having excessive synchronization betweenneurons is assumed to have a relationship with the small-worldarchitecture in the functional brain network. Ponten et al. conducted astudy doing graph analysis on intracerebral EEG recordings from patientshaving drug-resistant mesial temporal lobe epilepsy and found anincrease in the clustering coefficient in the lower frequency band (1-13Hz), and an increase in the path length in the alpha and theta bandsduring and after a seizure compared to interictal recordings. Thisimplies that the functional brain network changes from a more randomorganization to a small-world structure.

The efficiency of information transmission within a small-world networkimplies that some pathological damage to a brain may result in losingthe traits of a small-world network. In an Alzheimer's disease (AD)study conducted by de Hann et al., EEG data from patients with AD werecompared with control subjects. All subjects showed small-world networktraits in their functional graphs in all frequency bands, but in thebeta band alone, they observed significantly less small-worldness in theAD group compared to the controls.

Palva et al. did a study using magnetoencephalogram (MEG) and EEG toinvestigate the network properties when the subjects are engaged invisual working memory (VWM) tasks. The study recorded both MEG and EEGfrom 13 healthy subjects and applied PLV as the metric of associationbetween brain regions. Their results implied that small-world networkstructures appeared dynamically during the VWM task execution within thealpha and beta band. Moreover, the small-worldness was dependent on theVWM memory load.

Graph theoretic analysis has also been applied to the quantitative EEGdata of epilepsy patients. In a study by van Dellen, interictalelectrocorticography (EcoG) of 27 patients with pharmaco-resistanttemporal lobe epilepsy (TLE) were analyzed. Because an epileptic seizureis a manifestation of overly hyper-synchronous activity between neurons,a metric measuring synchronization, PLAI, was implemented on thoseelectrodes located on the temporal cortex and showed lower values onthose patients with a longer history of TLE. In addition, the clustercoefficient and small world index were both negatively correlated withTLE duration in the broad frequency band (0.5-48 Hz). This may haveresulted from the accumulative damage on the brain tissue due tointermittent excessive epileptic discharges over a long period of time.The optimal balanced small world structure had been impaired to form amore randomized one as the TLE duration continued.

In an interesting model undertaken by Rothkegel and Lehnertz, theyobserved the co-occurrence of local wave patterns and global collectivefiring in a two-dimensional small world network. Boersma et al.conducted a study on resting state EEG in developing young brains. Theyrecorded resting-state eyes-closed EEG (14 channels) from 227 childrenwhen they were 5 and 7 years of age and found out that the clusteringcoefficient increased in the alpha band with age. Path lengths increasedin all frequency bands with age. This suggests that a brain shifts fromrandom towards more ordered, small-world like configurations duringmaturation. Girls showed higher mean clustering coefficients in thealpha and beta bands compared with boys.

Schizophrenia has been suspected as the result of a more disconnectedbrain network among certain crucial areas in the brain. Rubinov et al.did a study investigating the “disconnection hypothesis.” They recordedresting state scalp EEG from 40 subjects with a recent first episode ofschizophrenia and another 40 healthy matched controls. Nonlinearinterdependences were identified from bipolar derivations of EEG dataand weighted graphs were generated. Graphs of both groups showedfeatures consistent with a small-world topology, but graphs in theschizophrenia group displayed lower clustering and shorter path lengths.This result can be interpreted as a pathological process that thesmall-world network transformed to a more randomized small-world networkin a schizophrenia brain. This randomization may be the reason whyschizophrenia evidences cognitive and behavioral disturbances. Inanother similar study, scalp EEG data from 20 schizophrenics and 20controls (age and sex matched) collected when they were performingworking memory tasks were analyzed using conventional coherence. Afterapplying a threshold on the values of coherence so that the number ofaverage edges is five in a graph, binary graphs were obtained. Theresults showed that schizophrenics bear a more random neuronal networkorganization, especially in the alpha band.

Inter-Hemispheric Power Asymmetry

The symmetric parts of a human body, such as limbs and sense organs,often have the same function and structure. The brain of a human alsohas symmetric shape although some functions happen on a dominant side.The EEG signals from a pair of symmetric channels usually have similarmorphologies. The asymmetry of EEG signals is regarded as a pathologicalconsequence or unusual phenomenon if the reason is not found. Manystudies have used asymmetry as an index to quantify the morbidity ofbrain disease or abnormal states. In this study, the hypothesis is thatthe relative frequency powers of symmetric pairs of EEG channels aremore different in CPS patients than that in PNES patients. The powers ofEEG signals in narrow frequency bands associate with different brainfunctions or motifs. The degree of asymmetry is quantified through therelative frequency power of several frequency bands and the T-index, astatistic measuring the degree of divergence between two groups.

Methods EEG Data Characteristics Subjects and EEG RecordingSpecifications

In this study, seven PNES and seven CPS patients were included. Allsubjects were 18 years of age or older who were admitted to the MedicalUniversity of South Carolina (MUSC, Charleston, S.C.) for inpatientseizure monitoring for diagnostic purposes or presurgical evaluation.The data collection procedure was approved by the Investigational ReviewBoards MUSC, and the Western Investigational Review Board (WIRB). TheEEG recordings at MUSC were obtained using XLTEK monitoring systems(Oakville, Ontario, Canada) with a sampling rate of 256 Hz and the EEGrecordings at AGH used 128-channel Nicolet BMSI-6000 systems (Viasys,Madison, Wis., USA) with a 400 Hz sampling rate. The EEGs recorded atboth institutions used a referential montage and the 19-electrodeinternational 10-20 system of electrode placement. The exact locationsof referential electrodes placed in the dataset were decided on-site andusually followed the recommended location of the Cz and Pz electrodes assuggested by the American Clinical Neurophysiology Society.

Awake and Relaxed State EEG Data Selection

The EEG signals of each subject were reviewed to select out the awakeand relaxed state sections containing the dominant alpha waves over theoccipital regions and no eye-blinking events around the frontalelectrodes. To classify patients using only interictal information, theawake and relaxed state EEG signals containing epileptiform dischargesor other suspicious epileptic activities were also excluded. Allselected awake and relaxed state sections were at least five hoursbefore the first CPS or PNES happened. During the awake and relaxedstate, the brain is supposed to be in a resting state and not activelyinvolved in any goal-oriented events. As a result, the awake and relaxedstates offer controlled background for the brain network to be comparedbetween subjects. On the other hand, an awake and alert section oflong-term continuous EEG from different patients could contain a varietyof ongoing psychological activity, raising concerns of confounding ifused for comparison analysis.

Functional Network Graph

Since oscillation frequency is a main characteristic of the brain andchanges when the brain undergoes different psychological conditions orexecutes different cognitive tasks, all selected epochs were filtered tospecific frequency bands including delta, theta, alpha, and beta. Allsignals were filtered using filters that do not distort the phaseinformation of the filtered signals so that the instantaneous phase timeseries of the original and filtered signals are the same. This iscrucial for this analysis because the connection strength was evaluatedusing phase information solely.

PLAI was estimated in every five-second epoch of all selected awake andrelaxed state sections. A weighted graph representing the functionalnetwork can be generated after all pair-wise PLAI are calculated amongstall electrodes. The weighted graph can also be presented as a weightedadjacency matrix as in FIG. 4-3. To convert the weighted graph to abinary one, a threshold can be chosen and applied to all PLAI valuessuch that the edges between vertices is either connected ordisconnected. The choice of threshold is done by controlling the densityof a graph so that the non trivial structure of the network can berevealed. It was selected to present the network in a graph with densityaround 0.75 so that only one quarter of the strongest connections(larger PLAI) remained in the functional network graph and the otherweaker edges (smaller PLAI) are ignored. The threshold in each epoch canbe different from each other and is increased slightly from a smallvalue until the desired density of a graph is reached. The thresholdchanges from epoch to epoch because the brain may undergo many phases ofsignal processing and show different connection strength betweenfunctional regions while the structure of the information flow shouldpersist in a small-world configuration so that the information isefficiently shared and processed among functional regions. For differentsubjects, it is reasonable to have individual thresholds for each epochso that the dominant structure of the functional network can be revealedand compared across subjects.

A network measure is a value quantifying a characteristic of thetopology of a network. Several measures have been proposed and used inanalysis of network structures. For determination of the existence of asmall-world configuration in a network, a clustering coefficient andminimum path length are two crucial measures to evaluate how small-worlda network is. A clustering coefficient quantifies how locally entangleda network is, and a minimum path length reflects how globally integrateda network is. These two measures of a functional network graph arecompared with those generated from 50 randomized graphs keeping the samenumber of vertices and edges. Given a clustering coefficient as Cp and aminimum path length as Lp for an epoch of selected awake and relaxedstate EEG; the randomized graphs generate Cp-s and Lp-s respectively.Given p_(i) is the number of edges between neighbors of vertex i, andk_(i) is the number of neighbors of vertex i, the local clusteringcoefficient can be calculated as Equation 4-14. In Equation 4-15, theclustering coefficient, Cp, is the average of local clusteringcoefficients of every vertices in the graph as V is the number ofvertices in the graph.

$\begin{matrix}{C_{i} = \frac{2p_{i}}{k_{i}\left( {k_{i} - 1} \right)}} & \left( {4\text{-}14} \right) \\{{Cp} = {\frac{1}{V}{\sum\limits_{i = 1}^{V}\; C_{i}}}} & \left( {4\text{-}15} \right) \\{{Lp} = {\frac{1}{V\left( {V - 1} \right)}{\sum\limits_{i,{j \in V},{i \neq j}}^{\;}\; m_{ij}}}} & \left( {4\text{-}16} \right)\end{matrix}$

In Equation 4-16, m_(ij) is the shortest path length from vertex i tovertex j. The ratio of Cp/<Cp-s> to Lp/<Lp-s> is called the small-worldnetwork index, λ. If λ>1, it is suffices to claim that the investigatedfunctional network is a small-world network because the investigatednetwork possesses a higher clustering coefficient or lower minimum pathlength comparing to 100 randomized networks possessing the same numberof edges and vertices.

Inter-Hemispheric Power Asymmetry

Frequency power density, X(f), can be estimated by applying discreteFourier transform on the interesting time series, x(t).

$\begin{matrix}{{X(f)} = {\sum\limits_{t = 1}^{N}\; {{x(t)}^{\frac{({{- 2}\; {\pi }})}{N}{({t - 1})}{({f - 1})}}}}} & \left( {4\text{-}17} \right)\end{matrix}$

However, discrete Fourier transform is not a consistent estimator andmodification is needed to better estimate frequency power density. Foreach five-second epoch, x(t) was divided into eight windows thatoverlapped half of the neighbor windows. For each window, hamming windowtechnique was applied to reduce the noise due to truncation of signaland power density was estimated. All power density functions from eightwindows were averaged as the estimate of the power density of the epoch.Relative powers of delta (1-4 Hz), theta (4-7 Hz), alpha (8-12 Hz), beta(13-30 Hz), and gamma (30-58 Hz) band were computed from the powerdensity function of frequency. For example, the alpha band relativepower would be the ratio of the sum of powers in the alpha band to thesum of powers from 1-58 Hz. These ratios, rp, were further transformedto a variable, rp′, as described in Equation 4-18 so that rp′ has adistribution close to the normal distribution

$\begin{matrix}{{rp}^{\prime} = {\log \left( \frac{rp}{1 - {rp}} \right)}} & \left( {4\text{-}18} \right)\end{matrix}$

The relative powers in each frequency band were compared to those of thecontralateral hemispheric EEG signals by the computation of the T-index.T-index is a function to compute the degree of divergence betweenpaired-samples from two groups. In this case, the sample groups were therelative powers of a certain frequency band from an anatomicallysymmetric, left and right, pair. Left and right relative powers in eachepoch should be paired up because they both reflected the state of thebrain during the same period of time.

$\begin{matrix}{{Tind}_{LR} = \frac{{\overset{\_}{D}}_{LR}}{{\hat{\sigma}}_{D_{LR}}/\sqrt{n}}} & \left( {4\text{-}19} \right)\end{matrix}$

In Equation 4-16, D _(LR) is the mean of the differences of relativepowers, rp′, of each awake and relaxed state EEG epoch and {circumflexover (σ)}_(D) _(LR) is the standard deviation of the differences. Due tothe fact that each subject has a different length of awake and relaxedstate EEG recordings, and the number of samples (degree of freedom), n,for calculating T-index should be fixed so that the comparison ismeaningful, random sampling was performed for each subject so that eachsubject has 36 (9 random samples from 4 continuous awake and relaxedstate EEG segments) epochs input to the T-index function.

Results Network Measures

For each five-second epoch, one functional network graph was generatedafter the PLAIs were estimated pair by pair. For each functional networkgraph, the clustering coefficient, minimum path length, and small-worldindex was calculated from the adjacency matrix. All network measures ofevery epoch from a patient were averaged into one value for each networkmeasure. As a result, every subject has one final value for each networkmeasure. Totally, the PNES or CPS patients group has seven subjects andtherefore, seven values for each network measure. The small-world index,λ, is larger than 1 in all frequency bands for both patient groups,supporting the existence of a small-world network structure in bothgroups in all frequency bands. A Student T-test was performed to test ifthe λ in CPS or PNES patient group was larger than one. The p-valuesshowed significance in all frequency bands for both groups (see Table4-1). Wilcoxon signed rank test was also performed and all subjectsshowed p-values as 0.0156. The Wilcoxon signed rank test results werethe same for all frequency bands since the sample number is always sevenand there were seven patients in both groups and all samples were largerthan one.

Following this, it was tested if these network measures showed adifference between the CPS and PNES patient groups. A Mann-Whitney Utest was performed to assess if the network measures were differentenough between both groups and the p-values were shown in the Table 4-2.Only the minimum path length ratio (Lp/Lp-s) in the delta band showed asignificant p-value.

Inter-Hemispheric Power Asymmetry

The inter-hemispheric power asymmetry was quantified by Tind_(LR). Inthe delta, theta, and alpha bands, the means of Tind_(LR) of every pair(besides T5-T6) in the CPS group were larger than that of the PNESgroup. In the beta band, the means of Tind_(LR) of every pair (besidesP3-P4) in the CPS group were larger than that of the PNES group. In thegamma band, the means of Tind_(LR) of every pair in the CPS group werelarger than that of the PNES group. Almost every anatomically symmetricpair in the CPS patient group showed more power asymmetry than that inthe PNES patient group. A Mann-Whitney U and Student's T test wereperformed to test if the means of divergences (T-indexes) ofindividually symmetric pairs were the same between the two groups andthe p-values are listed in Table 4-3 and Table 4-4 respectively. TheStudent's T test results showed that the relative powers of delta,theta, and alpha band on the C3-C4 pair in the CPS group hadsignificantly larger divergence than that in the PNES group.

Analysis of variance (ANOVA) partitions an observed variance intoseveral components of some possible factors and provides a test ofwhether the means of groups are equal. It was hypothesized that thevariance of the T-indexes, which quantify the inter-hemisphericasymmetry, can be partitioned into components explained by patientgroups, pairs and frequency bands. Two-way ANOVA considers two factorsin a linear model to explain the interesting dependent variable, andtests if the means of factor groups are the same. Amongst three factors,the main interest is to test if the patient group is a strong factorexplaining the variance of the T-index. The patient groups andanatomically symmetric pairs were first chosen as potential factors forexplaining the variance of T-indexes in the two-way ANOVA for eachfrequency band. The results showed that patient groups were thesignificant factor in the T-indexes, and the T-indexes weresignificantly different between the two patient groups under allfrequency bands except the beta band. The ANOVA p-values are presentedin Table 4-5. The patient groups and frequency bands were later used asfactors and did the two-way ANOVA again for each anatomically symmetricpair. The results are presented in Table 4-6 and show that the patientgroup was a significant factor explaining the variance of T-index forthose anatomically symmetric pairs in the frontal brain area.

Discussion

The global network structure and specific symmetric pair connectionswere investigated in both PNES and CPS EEG recordings. The small-worldnetwork index indicates how small-world a graph is comparing torelatively random graphs having the same number of vertices and edges.In Table 4-1, the small-world index was always bigger than 1 and theStudent's T-test showed significance in all frequency bands for bothpatient groups. The non-parametric Wilcoxon signed-rank test showed thesame results of being significant in all frequency bands. Both patientgroups showed small-world network structure in their functional graphsand the small-world indexes had no significant difference between twogroups. Only the minimum path length in the delta band showedsignificant difference between the two patient groups. Minimum pathlengths in the delta band of the CPS group were mostly shorter than thatof the PNES group (p-value=0.037). This implied that the globalintegration is more efficient in CPS patients than that in PNES patientsduring the awake and relaxed states. This is not surprising becauseseizure itself is a synchronous activity amongst many brain areas.Furthermore, about half of partial seizures happen during sleep.Moreover, temporal lobe complex partial seizures were more likely tosecondarily generalize during sleep than during wakefulness. My resultcould explain these interactions between partial seizures and sleep.Before entering the sleep state, the epileptic brain functional networkshowed higher global integration and could possibly facilitate thegeneration of seizures or the secondary generalization during sleep.This study did not include any sleep EEG data so the hypothesis shouldbe further pursued. Other network measures did not show that significantdifferences could result from several reasons and the possiblecombination of these reasons. First, the pathological network structuremay not manifest during the interictal awake and relaxed state of apatient. Second, the metric of association (PLAI) may not be sensitiveenough to differentiate the pathological nuance of neuronal interaction.Third, the network structure of the brain itself may be attack-tolerant.Even though there have been some lesions due to epileptic discharge, thebrain could be organized in a fashion such that tissue damage does notaffect much of the existence of the small-world feature in thefunctional network. These hypotheses can be confirmed through furtherstudy. For example, a study comparing the interictal and preictalnetwork structure could determine whether or not the pathologicalnetwork appears during the interictal state. Applying otherinterdependence measures could show a different perspective of thefunctional network in the awake and relaxed state. It would also beintriguing to design an integration index from different interdependencemeasures so that the integration index values correspond to thesynchronous activities of the greatest interest.

Inter-hemispheric power asymmetry is a more specific and local measureto investigate the quality of connection between hemispheres. For CPSpatients, the brain tissue around the foci could be damaged by therecurrent onsets of partial seizures. The lesion of the brain tissueshould affect the ensemble neuronal activity and cause the EEG signal todeviate from the EEG signal of the anatomically symmetric channel. Onthe other hand, PNES patients have more similarity between EEG signalsfrom anatomically symmetric channels due to symmetrically commensuratetissue integrity. Other than the lesion, the fixed focal hemisphere ofpartial seizures could result from substantially different pathologicalstructures within the hemisphere. If both of the abovementioned causesare valid and additive, the asymmetry could exacerbate.

TABLE 4-1 Small-world network index, λ, of functional networks of CPSand PNES patients during awake and relaxed state. Delta Theta Alpha BetaCPS Mean 1.182 1.119 1.124 1.108 SD 0.052 0.063 0.036 0.063 T-testp-value 8.72E−05* 0.002510* 0.000103* 0.003772* PNES Mean 1.150 1.1621.112 1.129 SD 0.054 0.049 0.038 0.037 T-test p-value 0.000315*0.000124* 0.000232* 9.96E−05* *significant result when the hypothesiswas that the λ is equal to one.

TABLE 4-2 Mann-Whitney U test results of network measures of CPS andPNES patient groups during awake and relaxed state. PNES CPS U test MeanSD Mean SD p-value λ Delta(1-4 Hz) 1.182 0.052 1.150 0.054 0.259Theta(4-7 Hz) 1.119 0.063 1.162 0.049 0.209 Alpha(8-12 Hz) 1.124 0.0361.112 0.038 0.456 Beta(13-30 Hz) 1.108 0.063 1.129 0.037 0.383 Lp/<Lp−s>Delta(1-4 Hz) 0.997 0.015 1.004 0.003  0.037* Theta(4-7 Hz) 1.003 0.0091.000 0.010 0.902 Alpha(8-12 Hz) 1.004 0.009 1.003 0.003 0.805Beta(13-30 Hz) 1.001 0.008 1.003 0.006 0.383 Cp/<Cp−s> Delta(1-4 Hz)1.177 0.041 1.155 0.056 0.383 Theta(4-7 Hz) 1.121 0.064 1.162 0.0430.383 Alpha(8-12 Hz) 1.128 0.035 1.115 0.040 0.535 Beta(13-30 Hz) 1.1090.062 1.132 0.039 0.318 *significant result.

TABLE 4-3 Mann-Whitney U test results of relative power asymmetry of CPSand PNES patient groups during awake and relaxed state. F3-F4 C3-C4P3-P4 O1-O2 F7-F8 T3-T4 T5-T6 Delta 0.209 0.017* 0.805 0.209 0.209 0.5350.053 Theta 0.097 0.073 0.805 0.535 0.165 0.038* 0.805 Alpha 0.318 0.0530.097 0.318 0.073 0.620 0.710 Beta 0.259 1.000 0.259 1.000 0.053 0.3830.902 Gamma 0.535 0.710 1.000 0.710 0.128 0.535 1.000 *significantresult.

TABLE 4-4 Student's two-sample T test results of relative powerasymmetry of CPS and PNES patient groups during awake and relaxed state.F3-F4 C3-C4 P3-P4 O1-O2 F7-F8 T3-T4 T5-T6 Delta 0.184 0.029* 0.293 0.2010.326 0.257 0.049* Theta 0.415 0.047* 0.633 0.412 0.236 0.054 0.944Alpha 0.262 0.050* 0.057 0.364 0.075 0.623 0.532 Beta 0.244 0.424 0.4500.644 0.110 0.170 0.752 Gamma 0.222 0.405 0.568 0.529 0.147 0.394 0.476*significant result.

TABLE 4-5 P-values of two-way ANOVA (patient groups as factor) in eachfrequency band. Frequency bands Delta Theta Alpha Beta Gamma P-value0.0036* 0.0034* 0.0019* 0.0584 0.0254* *significant result.

TABLE 4-6 P-values of two-way ANOVA (patient groups as factor) for eachanatomically symmetric pair. Symmetric pair F3-F4 C3-C4 P3-P4 O1-O2F7-F8 T3-T4 T5-T6 P-value 0.011* 0.001* 0.109 0.070 0.002* 0.017* 0.905*significant result.

Chapter 5 Connectivity Transition from Intericatl to Ictal States onNeocortical Epilepsy Preamble

Neocortical seizures originate in the neocortex—the external surfacepart of the cerebral hemispheres. Neocortical epilepsy differs frommesial temporal epilepsy in that it is difficult to clearly define asingle area from which the seizures originate. Since seizures associatedwith neocortical epilepsy generally do not respond well to medication,epilepsy surgery is often one of the few options that patients withneocortical epilepsy have. However, surgery for neocortical epilepsy hasa significantly lower success rate than other kinds of epilepsy. Inpatients with other types of epilepsy, such as mesial temporal lobeepilepsy (MTLE), surgeries can provide seizure freedom in more than 70to 90% of cases. Possible reasons that could result in the low successrate of neocortical epilepsy are: (1) there may be multiple seizureonset zones, (2) the initiation site may vary from seizure to seizure,and (3) with currently available technologies, it is difficult toprecisely identify the duration and extent of seizure onset zones. Arecent study by Lee et al. in 2005 on nonlesional neocortical epilepsyshowed that, based on the epileptogenic focus locations, only 47 out of89 patients were seizure-free (Engel Class I) after the surgery, and anadditional 7 experienced significant reduction in seizure frequency(Engel Class II). Therefore, developing a more reliable and effectivemethod for identifying suitable parts and critical regions forneocortical epilepsy surgery would be a major contribution to improvethe service of medicare and quality of life for those patients.

Many epilepsy researchers have postulated that there exists a network ofanatomical regions in the cerebral cortex and its connections thatinitiates the transitions between normal (interictal) and seizure(ictal) states, resulting in recurrent seizures. In this case,identifying the configuration of the network and understanding thecharacteristics and dynamics of the network could lead to a greaterunderstanding of neocortical epilepsy.

One of the most obvious hallmarks of the epileptogenic brain is thepresence of interictal spikes in the EEG. However, areas that generateinterictal spikes do not consistently correspond to the seizure onsetzones (brain site where the seizure discharge first appears).Nevertheless, it is very likely that areas with interictal spikes andthe seizure onset zones are both within the same epileptic brainnetwork. Furthermore, there may be other critical brain regions that arerelated to seizure occurrences but are not exhibited through interictalspikes or ictal discharges. My research, by analyzing intracranial EEGs,seeks to explore, quantitatively, the existence and structure of aneocortical epilepsy brain network. The methods and results sectionindicates how to identify the nodes and edges in the networkquantitatively, and how the characteristics of the brain network changefrom an interictal state to a seizure (ictal) state.

Another important diagnostic tool for neocortical epilepsy surgery ishigh-resolution magnetic resonance imaging (MRI). Often this offers ahigh predictive value for surgical outcomes. However, MRI does not workwell in 29% of patients with partial epilepsy. Many patients referred toepilepsy centers for surgery have normal MRI results (no lesion) andprevious studies report that surgical outcomes were poor for patientswith neocortical epilepsy with a normal MRI result. Intracranial EEGrecording is also indispensable for localizing seizure onset zones.However, sampling error can lead to false or missed localization duringintracranial EEG recording. This problem is particularly common in caseswith extratemporal seizure origin. Nevertheless, several studies havedemonstrated the usefulness of focus localization in neocorticalepilepsy with quantitative analysis on intracranial EEG. For example,Andrzejak et al. in 1999 applied methods derived from nonlinear dynamicson intracranial interictal EEGs (N=8) to lateralize the primaryepileptogenic area in two patients with neocortical epilepsy. They foundthat a nonlinear deterministic measure can contribute to an interictallocalization of the primary epileptogenic area in patients withneocortical epilepsy. Based on a “focus index” measure, Roopun et al.used a few examples to demonstrate that basic dynamic changes in focalepilepsy of neocortical origin may be useful in localizing the origin ofseizures. The findings of these studies and many others support thehypothesis that the dynamics of intracranial EEGs are associated withthe behavior of epileptogenic focus in neocortical epilepsy. Recentobservations in humans with MTLE and in the animal models for thiscondition helped conceptualize multifocal seizure onset (i.e. seizuresthat begin focally within different limbic structures with each seizure)as well as synchronized regional seizure onset (i.e. presumedsimultaneous seizure initiation). These observations, together withmultifocal physiological and anatomical changes in the animal modelshave raised the possibility of a widely distributed neural network(e.g., specific cortical and subcortical networks) in the gensis andexpression of partial onset seizures. Spencer (2002) defined a networkto be “a functionally and anatomically connected, bilaterallyrepresented, set of cortical and subcortical brain structures andregions in which activity in any one part affects activity in all theothers.” One central observation on which this definition is based isthat vulnerability to seizure activity in any one part of the network isinfluenced by activity everywhere else in the network, and that thenetwork as a whole is responsible for the clinical and electrographicphenomenon that may be associates with human seizures. She furtherdistinguishes between the network and the “seizure onset zone”. Theseizure onset zone for neocortical seizures is defined by the ictal EEG.It is the area where the seizure discharge is first seen on ECoG(subdural electrodes). It may involve one or many electrode sites. Insome patients, seizures (based on clinical manifestations) can startwith EEG discharges beginning in different “onset zones” and spreadingalong different routes as the seizure progresses.

Spencer (2002) suggested that the existence of an epileptic network issupported also by the responses to invasive therapy. If human epilepsyis the expression of specific, abnormally active, intrinsically definedand connected cortical/subcortical/bilateral networks, then one couldtheoretically alter seizure expression by intervening in any part of thespecific network. Operations involving anterior temporal lobe, medialstructures only, lateral structures only, or more or less extensivelateral temporal resection, can cure this disorder. Procedures with noanatomic overlap are similarly successful. This cannot be explainedunless the multiple areas are all critical in the production of theintractable seizures of this disorder. Then, interruption of the networkin any one of those areas would be (and apparently is) sufficient toalter the seizures. The similarly excellent response with cessation ofseizures after temporal lobe resection in well-selected patients whohave bilateral independent medial temporal lobe origin of seizures isanother example of the existence of a network, interference with whichat any site alters the expression of the intractable seizures.

Based on the above observations, it may be pertinent to consider studyof other kinds of phenomena in individual patients, which may define thenetwork in better terms than has been sought in the past because of thesingle-minded attention to defining regions of so-called seizure onset.For example, quantitative intracranial EEG analysis, backgroundpatterns, sleep effects on interictal and ictal activity, and othertypes of functional assessments may contribute considerably to theunderstanding of the role of networks in the expression of theepilepsies. Studying broad regions of brain structures related by thepresence of such networks, using quantitative EEG analyses andsophisticated approaches, may detect alterations in the behavior of thenetwork before the more traditional “seizure discharge” is seen, andallow prediction of seizures before manifestation clinically or ontraditional EEG.

Material

Intracranial EEGs provides the most convincing evidence to support thenetwork hypothesis. Because the entire network participates in theexpression of the seizure activity and can be entrained from any of itsvarious parts, initial electrical events (at “seizure onset”) may varyin their specific location of expression and occurrence within thenetwork. The initial area of apparent seizure involvement is not reallyan onset area, because “onset” could be expressed at any place in thenetwork, and might even vary from seizure to seizure in a given patient.This locational variability may produce different morphologies of“seizure onset” when EEG recording is performed in only one part of thenetwork. This may be the main reason for surgical failure in theneocortical epilepsy. A reliable quantitative method based onintracranial EEGs for determining the spatial distribution of the nodesof the epileptic network may lead to a better understanding of themechanisms that lead to the generation of a seizure, and provideinsights into more effective approaches to seizure control.

We selected one subject with neocortical epilepsy and used theintracranial EEG of the subject to construct a representative networkwhich is capable of showing the critical region with recordingelectrodes as vertices in the network. The patient is a 21-year oldmale, with a history of intractable seizures. Continuous VEEG monitoringwas performed with XLTEK monitoring systems (Oakville, Ontario, Canada)equipment with a sampling rate of 499.707 Hz. Totally, 44 electrodeswere applied on the left hemisphere cortex of subject. Data wereobtained, stored, and interpreted according to ACNS guidelines. Duringthe recording, frequent interictal spikes discharges were present at theG4, G5, G9, G14, G15, IH2, and IH3 electrodes.

Two EEG segments were cut and used to construct a brain network from acontinuous EEG recording containing one neocortical epilepsy seizurefrom the patient. The seizure events were identified by the recordingfacility based on both clinical observation and a review of the EEG. Inorder to avoid the potential overlapping between ictal, postictal andpreictal stages, only seizures that were at least 2 hours apart from theprevious or next one were used. The two segments analyzed were both 150min long. One contained a seizure starting at 120 min and the other wasseizure free. For the later, the previous or the next seizure happenedat least 5 hours away from the beginning or ending of the seizure freesegment. Both segments were preprocessed with a band-pass filter (1-220Hz) and several notch filters (60 Hz, 120 Hz, and 180 Hz) to remove thealias, power line artifacts and DC component.

Methods

This paper sought to identify the critical regions in the brain that areclosely related to the development of neocortical epileptic seizure. Inintracranial EEG recording, each electrode was treated (channel) as anode (vertex) in the network and cross-correlation functions were usedto calculate the weighting on the edge of each pair of nodes. To converta weighted network into a binary one, an edge stay was kept in thenetwork if and only if the two nodes adjacent to it were sufficientlyconnected.

The brain network was constructed as follows. Both 150-minute segmentswere further chopped into non-overlapping calculation windows of 6seconds to construct transient brain networks. In every calculationwindow, each electrode was treated as a node in the network. Pair-wirecross-correlation coefficients were calculated among all 44 nodes. Across-correlation coefficient, r(d), with delay d, was calculated asEquation 5-1.

$\begin{matrix}{{r_{x,y}(d)} = \frac{\sum\limits_{v = 1}^{N}\; \left\lbrack {\left( {{x(v)} - {mx}} \right) \times \left( {{y\left( {v - d} \right)} - {my}} \right)} \right\rbrack}{\sqrt{\sum\limits_{v = 1}^{N}\; \left( {{x(v)} - {mx}} \right)^{2}}\sqrt{\sum\limits_{v = 1}^{N}\; \left( {{y(v)} - {my}} \right)^{2}}}} & \left( {5\text{-}1} \right) \\{{{mx} = {\frac{1}{N}{\sum\limits_{v = 1}^{N}\; {x(v)}}}},{{my} = {\frac{1}{N}{\sum\limits_{v = 1}^{N}\; {y(v)}}}}} & \left( {5\text{-}2} \right)\end{matrix}$

where x(v) and y(v) are EEG signal values of a pair of electrodes attime v and N is the number of sampling points in each 6-secondcalculation window. The cross-correlation coefficients with variousdelay d (range from −1 to +1 second) were calculated. Then, the maximumof the absolute value of the cross-correlation coefficient was obtainedby selecting delay d. For each pair of nodes p and q, the weight of edgee_(p,q) is defined by this maximum value. In this way, one obtains aweighted network. To convert it to a binary network, only the edges witha weight higher than a threshold (0.9) stay in the network.

$\begin{matrix}{e_{p,q} = {\max\limits_{d}\left\lbrack {r_{p,q}(d)} \right\rbrack}} & \left( {5\text{-}3} \right) \\{E_{p,q} = {H\left( {e_{p,q} - {thr}} \right)}} & \left( {5\text{-}4} \right)\end{matrix}$

In each six-second epoch, a graph representing the instant brain networkwas constructed and then the degree of electrode x, D_(x), wascalculated in Equation 5-5.

$\begin{matrix}{D_{x} = {\sum\limits_{y = 1}^{N_{nodes}}\; E_{p,q}}} & \left( {5\text{-}5} \right)\end{matrix}$

In Equation 5-5, N_(nodes) is the total number of nodes implanted underthe scalp and equals 44 for this neocortical epilepsy patient. Thefive-minute mean occurrences of degree of node x, D_(x(5-min)), wascalculated and then an average degree over all nodes, AveD wascalculated by averaging D_(x(5-min)) over all channels as Equation 5-6.

$\begin{matrix}{{AveD} = {\frac{1}{N_{nodes}} \times \left\lbrack {\sum\limits_{x = 1}^{N_{nodes}}\; {D_{x}\left( {5 - \min} \right)}} \right\rbrack}} & \left( {5\text{-}6} \right)\end{matrix}$

Results and Hypothesis Testing

After the coupling strength between each pair of nodes was quantified asthe maximum cross-correlation, the functional epileptic network wasconstructed. There are many features of network structures, and degreeis one of the most widely-used basic properties of network analysis.Based on the degree results, the connectivity of an instant network wasdefined as the average degree of all nodes during a five minute periodand quantified as AveD. One of the main purposes of this study was tofind if there was a significant change in the instant network structuresbefore and after a seizure onset. Therefore, one of those properhypotheses to test this idea was that the connectivity of an epilepticbrain network increases when a brain approaches a seizure onset and thendecreases after the onset. The results are shown as follows and thehypothesis test will be done in the discussion section.

To have a better time resolution and globally spatial layout ofconnectivity, AveD was used to present the connectivity. After the AveDcalculation, the trajectories of AveD were drawn. The slope of theaverage degree over all channels was estimated between the beginning and120 minute time point (i.e., onset time) of the segment. The slope fromthe 0 to 120 minute time point in the seizure segment was 0.02/5 minutes(p=0.028), and was −0.0036 (p=0.447) during the same time period of theinterictal segment. The degrees of nodes in the network were averagedover 300 calculation windows that spanned over 30 minutes, and thenplotted in FIG. 5-4.

Discussion

In FIG. 5-2 and FIG. 5-3, the regression slope of the segment with anonset is larger than zero, while the segment without any onset is lessthan zero. This suggested that the connectivity increases before anonset. Although the interictal segment had high values of AveD aroundthe 95 minute, the regression slope was not larger than zero. Based onthis observation, the following can be concluded. If the recording isdone during interictal state, the connectivity should fluctuate within arange and have a regression slope around zero. If one has a preictalrecording, one should more likely detect a positive regression slopewithin a period right before the seizure onset. The p-value providesinformation about how unlikely the regression slope is compared to thevalue of zero. The regression slope equal to zero means that theconnectivity does not correlate with the time within the segment. Thisis also what was presumed to happen during an interictal state.Additional support for this hypothesis is the result that the p-value isequal to 0.0028 (<0.05) in the segment preceding an onset and 0.447(>0.05) in the interictal segment.

An epileptic seizure can be interpreted as an activity during which theneurons from different region overly correlate with each other.Therefore, one should observe a decrease of connectivity after aseizure. Network connectivity in terms of AveD 30 minutes before and 30minutes after the seizure onset were compared. The average degree beforethe seizure was 0.33 versus 0.09 after the seizure (p=0.027). Applyingthe same comparison of average degrees in the interictal segment, thedifference was not significant (p=0.366). The results of this pilotstudy suggest that transitions in the brain networks may exist and arerelated to the underlying dynamics of seizures caused by neocorticalepilepsy. An increase of connectivity can be observed before seizureonset. The onset of a seizure can cause the reset of the connectivitycloser to a baseline level.

Prokopyev et al. (2003) conducted additional analyses on the brain andfound significant changes in the network characteristic during aseizure. In this paper, a decrease of connections was observed inresults after the onset (in FIG. 5-4); this is consistent with the resetphenomenon after seizure. Monto et al. in 2007 studied epileptic brainnetworks by quantifying the long-range temporal correlation in subduralhuman EEG recordings. Their observations on the spread of abnormallylarge LRTCs suggested that the epileptic focus is associated withsignificant changes in network behavior even in the cortical areasimmediately surrounding the clinically determined focus. By applyingsynchronization and graph analysis to intracranial EEG recordings,Ponten et al. (2007) investigated the hypothesis that functionalneuronal networks during temporal lobe seizures change in configurationbefore and during seizures. The study found that the functional brainnetworks change from a more random configuration during an interictalrecording period to a more ordered configuration during seizures,especially during seizure spreading. The authors further suggested thatthe findings supported the theory that a random network (duringinterictal periods) even had a stronger tendency to synchronize, whichcould cause seizures. More recently, Zaveri et al. (2009) calculated amagnitude squared coherence on background intracranial EEG as a measureof functional connectivity to investigate the network effects within andoutside the seizure onset area. The analysis demonstrated an inverserelationship between the connectivity strength and the distance from theseizure onset area (especially during the frequency band). Anothersimilar work by Ortega et al. in 2008 showed the presence of clusters ofincreased synchronization in different locations on the lateral temporalcortex in patients with temporal lobe epilepsy. As yet unproven, it ispossible that increased connectivity (as can be defined by coherence,synchronization clusters, or nonzero functional connectivity) helpsinitiate seizures. If increased brain connectivity or alteration innetwork topology helps initiate seizures, then network nodes andpathways could serve as targets for resective or disconnective surgery,implantable devices, and investigations of seizure anticipation.Applying classification techniques for circumscribing some communitieswithin the network could be done after the network is constructed andthese communities may serve as good entities for observing theentrainment preceding a seizure. This is due, in part, because they aresupposed to be less correlated during interictal periods and thenentrain as evolving to ictal state. Several methods of identifyingcommunities in network were introduced in Porter's work in 2009.

Chapter 6 Conclusion

My work investigated several time series features applied to severalapplications relating to epileptic EEG signals. The first part of myresearch focused on a possible scheme for seizure prediction. Whilethere have been numerous studies that show promising results inanticipating seizure onsets using intracranial EEG, transferring thesesuccesses to seizure prediction based on scalp EEG is still achallenging task. Many studies so far have investigated the seizureprediction problem only under very well-controlled conditions.Developing a robust scheme for seizure prediction outside the laboratoryenvironment requires a variety of signal processing techniques. Thefirst part of my research explored the performance of a singlefeature-based algorithm on long-term continuous scalp EEG. The focushere was on evaluating an automated seizure prediction algorithm thatissues seizure warnings by monitoring the convergence of signalregularity among EEG channels in continuous long-term scalp EEGrecordings. For each trial, the algorithm was optimized with a trainingdataset and its' performance was evaluated using a separate testdataset. In the test dataset, the algorithm achieved an averagesensitivity of 65% (i.e., ˜2 out of every 3 seizures) with an averagefalse positive rate of 0.24 per hour (˜1 per 4 hours). The algorithmperformance in the test dataset was nearly the same as that in thetraining dataset. This implies that the algorithm generated a stableperformance in epileptic patients. Statistically, the overallsensitivity achieved by the test algorithm was better than a randomwarning scheme (p-value=0.0145). While these results are encouraging,the performance of the algorithm may not yet be sufficient for someclinical applications. For example, the performance would be of limitedutility for inpatient monitoring applications, such as the EMU setting,due to the high false positive rate. However, the performance may beuseful for driving seizure control devices, such as the vagus nervestimulator. Furthermore, in this type of application, fine tuning theprediction algorithm to optimize performance for each individual patientcould yield even better performance. Although the performance was notyet sufficient for its usage in a clinical environment, the resultsdemonstrated the supremacy of the algorithm compared to a random scheme.This implied that the method found some events preceding the onsets ofseizures. Further analysis needs to follow to be more specific about thepositive findings in this study.

It is also informative to investigate the common characteristics ofthese negative results. For instance, the false-positives may resultmostly from a certain type of event and could be muted if the triggeringevent is found. This could bring down the false-positive rate andpossibly further increase the sensitivity of the algorithm when theparameters are chosen to be more sensitive and the false-positive rateis well suppressed to an acceptable level. Incorporating engineeringtechniques could help the prediction scheme become more applicable inreality. For a complex system such as the brain, it may be too ambitiousto hope that only one dynamic feature or method would be sufficient toaccomplish effective seizure prediction. Aggregating the advantages fromseveral time series dynamic features can possibly join the benefit ofeach feature. Until now, most studies have compared one kind of methodto another, such as linear vs. nonlinear and univariate vs. bivariatefeatures. These results demonstrated that the bivariate featuresoutperformed univarite features. Nonlinear features have been popularsince 1999 but progress since then has not been substantial. Untilrecently, linear methods such as the AR model or coherence have stillbeen used and have generated satisfactory results in many EEG researchfields including seizure prediction. Different features explain the samesignal from different perspectives and should be clearly correlated tothe raw signal so that one can better utilize a feature at the propertime for a suitable purpose. Besides signal features, additionalengineering methods could benefit the prediction scheme in thepreparation and post-feature stages. In the preparation stage,filtering, inverse source locating, and network analysis could helpextract more seizure-relevant information from the EEG signals beforethe dynamic features are applied to the signals. During the post-featurestages, a decision to issue a warning or an index denoting thesusceptibility of having a seizure should be output. A decider should bedesigned to consider situational information about the differentfeatures from the patient and integrate them in an optimal way.Optimization skills could be involved to assist in the decision process.Additionally, a prediction scheme with an updating threshold or alearning classifier would be more adaptive to the possibly changingbackground of EEG activities if one considers applying the predictionscheme in an ambulant environment. There are still many aspects that canbe improved and deserve greater analysis or hypothesis testing inseizure prediction. The understanding of applicable features on the rawdata and the integration of information are fundamental to accomplishthis challenging task.

In my initial study, the focus was on the temporal transition from theinterictal to ictal states of a group of epileptic patients. Followingthis, the possible differences of interictal network features betweenpsychologically and physiologically epileptic patient groups wasinvestigated. Graph theoretical analysis was applied to reveal thefunctional network structure of CPS and PNES patients when they were inthe awake and relaxed interictal state. The results showed the existenceof a small-world network configuration in both patient groups.Furthermore, the delta-band minimum path length was significantlysmaller in the CPS patient group. This implied that the low frequencyglobal integration was more efficient in CPS than in PNES patients. Bothgroups had the same small-world network indexes and clusteringcoefficients. These results implied that the brain either somehowmaintained the structure of a small-world network or the originalsmall-world network was attack-tolerant so that the pathologicalinfluence of epileptic discharges did not bring down the small-worldconfiguration. The results above viewed the connection in acomprehensive scope over the whole brain. However, if the pathologicaleffect was only locally misleading, the large-scale analysis could failto track down the difference. To take a closer look at the connectionstrength of several specific pairs, the frequency power information wasretrieved and compared the inter-hemispheric asymmetry for both groups.The difference of power asymmetry was revealed in all frequency bands onmost of the frontal pairs. The results showed that CPS patients had morepower discrepancies between contralateral pairs than PNES patients.These differences were not significant for those pairs around the rearpart of the brain which may result from the epoch selection criterion.The awake and relaxed state EEG segments were selected based on thedominant alpha-band oscillations appearing around the occipitalchannels. This selection criterion forced the presence of high alphapower in all segments and may have reduced the power of the statisticaltest. Some EEG studies in other fields also focus on power asymmetryonly for the frontal region. Future research for this study couldinvolve applying a classifier to distinguish if a subject is a CPS orPNES patient by analyzing only a short period of the awake and relaxedstate EEG. This would greatly decrease the recording time and resourcesof an EMU. Due to the requirement that a CPS patient should receiveanti-epileptic treatment, the classifier should be designed to precludemisdiagnosing a CPS patient as a PNES patient, but could be permitted toincorrectly identify a few PNES patients as CPS patients. The results ofthis study also suggest that future hypotheses may need to include moresensitive scales for detecting these important distinctions. For the bigscale network analysis, these differences were not obvious orsignificant. However, a more precise comparison of specific pairsdisclosed more locally detailed information and showed the expecteddifference as hypothesized.

Therefore, at least the following is claimed:
 1. A method for cerebraldiagnosis, comprising: obtaining a plurality of electroencephalogram(EEG) signals from a plurality of sensors positioned about a scalp of asubject; conditioning data from the plurality of EEG signals to removeartifacts; generating a cerebral network model based at least in partupon the conditioned EEG signal data; determining network features basedupon the cerebral network model; and determining a cerebral condition ofthe subject based at least in part upon the network features.
 2. Themethod of claim 1, further comprising determining EEG signal featuresfrom the conditioned EEG signal data, wherein the cerebral condition isbase at least in part upon the EEG signal features and the networkfeatures.
 3. The method of claim 2, wherein the cerebral network modelis generated based at least in part upon the determined EEG signalfeatures.
 4. The method of claim 1, wherein the removed artifactsinclude eye movement artifacts and muscle movement artifacts.
 5. Themethod of claim 1, wherein the removed artifacts include sensor relatedartifacts.
 6. The method of claim 1, wherein generating the cerebralnetwork model comprises: generating a weighted graph based upon EEGsignal features determined from the conditioned EEG signal data; andconverting the weighted graph to a binary graph based upon a predefinedthreshold.
 7. The method of claim 1, wherein determine network featurescomprises determining global network characteristics.
 8. The method ofclaim 1, wherein determine network features comprises identifying hubsof the cerebral network model.
 9. The method of claim 1, whereindetermining the cerebral condition of the subject comprises determininga location of an abnormal condition.
 10. The method of claim 1, whereindetermining the cerebral condition of the subject comprises determininga severity of an abnormal condition.
 11. A method for cerebraldiagnosis, comprising: positioning a plurality of electroencephalogram(EEG) sensors about a scalp of a subject: determining a recordingcondition for each of the plurality of EEG sensors based upon predefinedsensor criteria; in response to an unacceptable recording condition forat least one of the plurality of EEG sensors based upon the predefinedsensor criteria, providing an indication of the EEG sensor correspondingto the unacceptable recording condition and provide procedures tocorrect the unacceptable recording condition.
 12. The method of claim11, further comprising in response to acceptable recording conditionsfor the plurality of EEG sensors based upon the predefined sensorcriteria, obtaining a plurality of EEG signals from the plurality of EEGsensors.
 13. The method of claim 12, further comprising amplificationand filtering of the obtained plurality of EEG signals.
 14. The methodof claim 12, further comprising sampling the obtained plurality of EEGsignals to obtain digital EEG data.
 15. The method of claim 14, whereinthe digital EEG data is stored in memory.
 16. A system for cerebraldiagnosis, comprising: an electroencephalogram (EEG) recording moduleconfigured to acquire signals from a plurality of sensors positionedabout a scalp of a subject; a signal conditioning module configured tocondition EEG signal data from the plurality of EEG signals; a signalanalysis module configured to determine EEG signal features and cerebralnetwork features based at least in part upon the conditioned EEG signaldata; and a condition classification module configured to determine acerebral condition of the subject based at least in part upon thedetermined features.
 17. The system of claim 16, further comprising anelectrode application module configured to verify a recording conditionof each of the plurality of sensors based upon predefined sensorcriteria.
 18. The system of claim 16, wherein the signal conditioningmodule is configured to remove artifacts associated with movement of thesubject from the EEG signal data.
 19. The system of claim 16, whereinthe signal analysis module is configured to generate generating acerebral network model based at least in part upon the conditioned EEGsignal data.
 20. The system of claim 16, wherein the conditionclassification module is configured to identify location and severity ofan abnormal condition.